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ABSTRACT 

We present traphic, a novel radiative transfer scheme for Smoothed Particle Hydro- 
dynamics (SPH) simulations, traphic is designed for use in simulations exhibiting a 
wide dynamic range in physical length scales and containing a large number of light 
sources. It is adaptive both in space and in angle and can be employed for applica- 
tion on distributed memory machines. The commonly encountered computationally 
expensive scaling with the number of light sources in the simulation is avoided by in- 
troducing a source merging procedure. The (time-dependent) radiative transfer equa- 
tion is solved by tracing individual photon packets in an explicitly photon-conserving 
manner directly on the unstructured grid traced out by the set of SPH particles. To 
accomplish directed transport of radiation despite the irregular spatial distribution 
of the SPH particles, photons are guided inside cones. We present and test a parallel 
numerical implementation of TRAPHIC in the SPH code GADGET-2, specified for the 
transport of mono-chromatic hydrogen-ionizing radiation. The results of the tests are 
in excellent agreement with both analytic solutions and results obtained with other 
state-of-the-art radiative transfer codes. 

Key words: methods: numerical - radiative transfer - hydrodynamics - cosmology: 
large-scale structure of the Universe - HII regions - diffuse radiation 



1 INTRODUCTION 

Radiation is one of the fundamental constituents of our Uni- 
verse. Its interaction with baryons may lead to an energy 
exchange that can both heat and cool the matter, initiating 
pressure forces that may strongly influence the subsequent 
hydrodynamical evolution. Radiation may also exert a di- 
rect pressure force upon the matter through the exchange 
of momentum. Radiative interactions are furthermore of- 
ten the dominating process in governing the excitation and 
ionization state of atoms and molecules. The inclusion of 
the transport of radiation into hydrodynamical simulations 
may therefore provide the key for interpreting the outcome 
of physical experiments and observational campaigns. 

To perform hydrodynamical simulations, the La- 

grangi an technique Smoothed Particl e Hydrodynamics 

(SPH; iGingold fc Monagnanlll977l : Illucvlll9771 ) is often em- 
ployed. In SPH, the continuum fluid is discretized using a 
finite set of point particles, each carrying its own collection 
of variables. The deformation of the fluid by internal and 
external processes manifests itself in a steady redistribution 
of the point particles in space. All it takes to determine the 
values of physical field variables at a given point in the sim- 
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ulation box is to perform a weighted average over the values 
the relevant variables take on the particles in its surround- 
ing. This elegant simplicity of the SPH technique is one of 
the many reasons for its success. 

Although the first radiative transfer calculation was al- 
ready included in SPH at the very bir th of this n umerical 
technique more than thirty years ago (|Lucy|[l977h . the de- 
tailed treatment of radiation transport in SPH is still an 
enormous challenge. One of the main reasons for this is cer- 
tainly the high dimensionality of the problem. In fact, the 
radiative transfer equation depends on no less than seven 
variables (three space coordinates, two angles, frequency and 
time). Moreover, existing numerical schemes to solve the ra- 
diative transfer equation have most often only been formu- 
lated for use with uniform grids. Despite these difficulties, 
there has been encouraging progress over the last years, with 
several interesting and rather different approaches (see Sec- 
tion [2]). In this paper we present a new method to solve 
the radiative transfer equation in SPH. We specifically de- 
signed our method for use in hydrodynamical simulations 
exhibiting a wide dynamic range in physical length scales 
and containing a large number of light sources. 

A prominent example for such simulations are cosmo- 
logical simulations of large-scale structure formation. Per- 
forming cosmological simulations is an exceedingly demand- 
ing computational task. Difficulties in describing the growth 



2 Andreas H. Pawlik & Joop Schaye 



of the initially tiny matter density perturbations produced 
before and during the event of recombination and their 
metamorphosis into the rich structure observable today do 
not only arise because of our ignorance of how to properly 
model the governing physical processes. Often, it is simply 
the lack of computational power which prevents us from 
faithfully representing the basic actions involved: Cosmolog- 
ical simulations are both time consuming and memory ex- 
haustive. To overcome these computational challenges, one 
can make use of advanced techniques and resort to clever 
approximations, reducing the computational effort. 

Consider the wide range of scales encountered in cos- 
mological hydrodynamical simulations. According to the hi- 
erarchical model of structure formation, the first structures 
and building blocks of the evolving universe are expected 
to form at small scales. The non-linear evolution at these 
scales shapes the distribution of the matter at all times, 
thus necessitating simulations of high spatial resolution. On 
the other hand, we also require sufficiently large simulation 
boxes in order to properly account for the modulation of the 
small-scale nonlinear evo lution by the large-sca le structure 
formation processes (e.g. Ba rkana fc Loebll2004 ) and not to 
be deceived by the cosmic scatter. 

To accommodate these two antithetic demands while 
keeping the number of particles representing the matter low 
enough to be computationally manageable, spatially adap- 
tive SPH simulations have been invoked. It is then imme- 
diately clear that when solving the transport of radiation 
along with the gravito-hydrodynamics of the matter, one 
requires the radiative transfer scheme to be adaptive, too. 
Even spatially adaptive cosmological SPH simulations, how- 
ever, make use of hundreds of millions of particles and are 
therefore still extremely memory-consuming. It is then indis- 
pensable to distribute the computational load over a large 
number of machines. For this reason, we require a radia- 
tive transfer scheme that is parallel on distributed memory 
machines. 

When performing radiative transfer simulations, 
sources of light are assigned a special importance. As a re- 
sult, the computation time of most of the available radiative 
transfer schemes scales linearly with the number of sources 
in the simulation box. However, in cosmological simulations, 
even at times as early as 1 billion years after the Big Bang, 
i.e. at a redshift of z ~ 6, a non-negligible amount (~ 1 per 
cent) of baryonic matter has undergone star formation. In 
addition to these stellar sources of light, the intergalactic 
gas emits photons too, producing a radiation component of- 
ten referred to as diffuse radiation. Hence, without breaking 
the linear scaling with the number of light sources, radiative 
transfer simulations at the resolution of state-of-the-art cos- 
mological hydrodynamic simulations need to be dispatched 
to the realm of the future. 

In this paper we present a radiative transfer scheme 
for use in SPH simulations that is adaptive, parallel on 
distributed memory and that avoids the linear scaling of 
the computation time with the number of sources, mak- 
ing it ideal for application in large simulations covering a 
wide range of length scales and containing many sources. 
In our scheme we follow the propagation of individual 
photon packets. Hence, it is explicitly photon- conserving 
(|Abel. Norman, fc Madaulll999T l and can be applied to solve 
the time- dependent radiative transfer equation. Because the 



photon packets are traced in cones, we refer to our scheme 
as traphic - TRAnsport of PHotons In Cones. The intro- 
duction of cones is required in order to perform the trans- 
port of photon packets directly on the unstructured grid 
defined by the SPH particles. Although we have designed 
our radiative transfer scheme to be readily coupled to the 
hydrodynamic evolution of the matter in SPH simulations, 
here we limit ourselves to the description and testing of the 
radiative transfer scheme itself. We will report on fully cou- 
pled radiation-hydrodynamics SPH simulations employing 
traphic in future work. 

The outline for the rest of this paper is as follows. We 
start with a brief review of the principles of SPH that we 
consider crucial for the understanding of the present work. In 
Section [3] we then recall the radiative transfer equation and 
place our radiative transfer method in context by reviewing 
the approaches that have been used to solve the radiative 
transfer problem in SPH so far. In Section U we give a de- 
tailed description of the ideas behind our method. Through- 
out we will emphasize how we satisfy the requirements set 
by our primary aim of performing radiative transfer in sim- 
ulations exhibiting a wide range in length scales and a large 
number of light sources. Thereafter, in Section [5] we apply 
our method to the transport of hydrogen-ionizing radiation, 
describe its numerical implementation in the parallel state- 
of-the-art SPH code GADGET- 2 and report its behaviour in 
test problems. Finally, we summarise our approach and con- 
clude with an outlook. 



2 SMOOTHED PARTICLE HYDRODYNAMICS 

Excellent reviews of SPH e xist (see e.g. iMonaghanl 120051 ; 
lMonaghan|[l992l ; iPricd 120051 ). Here, we just briefly outline 
the basic principles we consider critical for the understand- 
ing of our radiative transfer method. 

At the basis of SPH lies the representation of a field 
A(r) by its integral interpolant Ai(r), 



Ai(r) = / dV A(T')W(r-r',h), 



(1) 



where the smoothing length h determines the spatial reso- 
lution. The interpolation kernel W satisfies 



j dV W(r - r', h) 
lim W(r - r', h) 

h->0 



5 D (v- 



(2) 
(3) 



where 8d is the Dirac delta function, such that Ai(r) co- 
incides with A(r) in the limit h — > 0. The last two condi- 
tions do not fix the functional form for W. An often adopted 
choice is the spherically symmetric compact spline 



W{T-v',h) 
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(4) 



where r — |r — r'|. From here on, when referring to the 
interpolation kernel W , we assume that it is of this form. 
The discrete representation of a fluid by SPH particles is 
achieved by approximating the integral interpolant (Eq. [1} 
by the summation interpolant, 



Ai(r) « A s (r) 



A{Tj)W{T-Ti,h), 



(5) 
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where the summation is over the particles of mass m and 
mass density p. For a self-contained numerical treatment, 
any field needs to be discretized only at the positions of the 
particles i representing the fluid, 



Ai = A s (n) =r^A(r 3 )%- r it A). 

A Pi 



(6) 



Since the kernel W is compact, only local information needs 
to be accessed for the evaluation of the last sum, which 
can be carried out in a computationally efficient manner 
in parallel on distributed memory machines. 

For the interpretation of Eq.|6] two m ain approaches can 
be ta ken, depending on the choice for h l|Hernquist fc Kata 
1989). In the scattering approach each particle j is consid- 
ered as a cloud of radius hj and contributes to the field at 
the position of particle i with the weight W(ri — ij, hj). In 
the gathering approach, each particle i searches the sphere 
of radius hi (in the following referred to as the sphere of in- 
fluence, or neighbourhood) for particles (its neighbours) and 
obtains an estimate of the field value at its position by sum- 
ming the field values at the positions of the neighbours j, 
each weighted by W{vi — Tj,hi). 

The two interpretations are identical only if the spatial 
resolution is fixed throughout the fluid, i.e. if hi = hj = h. 
However, the SPH formalism only unfolds its true strength 
by allowing the resolution to vary in space, according to the 
local density field. This is usually achieved by either directly 
fixing the number of neighbours N ng b for all particles, or by 
requiring that the kernel vo lume contains a constant mass 
jSpringel fc Hernquisdf2002h . 



3 RADIATIVE TRANSFER IN SPH - 
PREVIOUS WORK 

Before describing traphic, our method to solve the radia- 
tive transfer equation in SPH, we briefly recall the radiative 
transfer equation and give a short overview of the main nu- 
merical methods that have been employed so far to obtain 
its solution in SPH simulations. 

The classical equation of radiat ive transfer reads (see 
e.g. iMihalas fc Weibel Mihala3ll984h 

' <>r " + n ■ VJ„ = t v — Kvplv (7) 



c dt 
In Eq. /, 



I v (v, n, t) is the monochromatic inten- 



sity (units ergs cm s 



Hz" 



_1 ) of frequency v 



at position r, n is a unit vector along the direction of 
light propagation and c is the speed of light. Sources and 
sinks of radiation are described by the emissivity e„ (units 
ergs cm -3 s _1 Hz -1 sr _1 ) and the mass absorption coeffi- 
cient K u (units cm 2 g _1 ), respectively. In general, both e v 
and n v are functions of r, n and time t. If we define the pho- 
ton number density ip u such that ipudfldv is the number of 
photons per unit volume with frequencies (y, v + dv) trav- 
elling with velocity c into a solid angle dQ around n, the 
intensity is given by I v = ch v vi\) v , where h p is the Planck 
constant. Eq. [7] is a partial differential equation depending 
on seven variables, which in general is hard to solve. Differ- 
ent approximations have been employed to obtain its solu- 
tion, giving rise to different numerical approaches. Below we 
discuss the main approaches that have been taken so far to 
accomplish the transport of radiation in SPH simulations. 



In his study of optically thick proto-stars, iLucvl (fl977h 
modelled the transport of radiation as a diffusion process 
by including a corresponding term in th e SPH formula- 
tion of the energy equation. iBrookshawl \ 19941 ) p ointed 
out d rawbacks of the particular formulation used in ILucvl 
(1977), and re-phrased the diffusion equation in SPH such 
as to reduce its sensitivity to particle disorder. This new 
formulation of the d i ffusion equation was employed by 
IViau. Bastien. fc Chal (2006) to perform collapse simula- 
tions of centrally condensed clouds. 

Treating the transport of radiation in the diffusion limit 
is, however, only a good approximation in optically thick 
media. In the optically thin regime, infinite signal propaga- 
tion speeds result, since the diffusion equation is a partial 
differential equation of parabolic type. This defect can be 
remedied by supplementing the diffusion equation with a 
so-called flux-limiter. SPH simulations with flux-lim it ed di f- 
fusion have been ca rr ied o ut by e.g. iHerant et al.l lll994h. 
IWhitehouse fe Bate! (fioop^. iFrver. Rockefeller, fc Warrenl 



1 20061 ) and Mayer et al.l 1 20071 ') . covering a wide range of 
physical applications, from supernovae explosions and proto- 
stellar collapse to the study of proto-planetary disks. The 
(flux-limited) diffusion approach to solve the radiative trans- 
fer equation fails in complex geometries. In particular, 
opaque obstacles illuminated by a point source do not cast 
sharp shadows, because the shadow region is filled in by the 
diffusion. 

To solve the radiative transfer equation without be- 
ing restricted to the optically thick regime or simple 
geometries, Monte Carl o techniques can be employed 
JOxlev fc Wooffsonl Eool. IStamatellos fc Whitworthl 120051 . 
ICroft fc Altavl 120071 . ISemelin. Combes, fc Baekl \200H) . In 
Monte Carlo radiative transfer, individual photon packets 
emitted by each source are directly followed as they pass 
through the matter, thus simulating the physical process of 
radiation transport as encountered in nature. While Monte 
Carlo simulations allow realistic shadows to be cast behind 
opaque obstacles, they are computationally very demanding, 
since the Poisson noise inherent to the statistical description 
of the radiation field leads to a signal-to-noise ratio that 
grows only with the square-root of the number of photon 
packets emitted. 

Ray-tracing schemes keep the advantages of Monte 
Carlo radiative transfer simulations, whilst not being af- 
fected by the statistical noise. In short, in ray-tracing 
schemes rays are cast from each source throughout the sim- 
ulation box. The optical depth to each point in space is 
calculated along theses rays, and the attenuation of the flux 
emitted by the sources is obtained. Variants of ray-tracing 
working with photons instead of fluxes blur the differences 
with Monte-Carlo schemes. 

Ray-tracing schemes are most easily implemented 
on a regular grid superimposed on the SPH den- 
sity field. However, this implies that all information 
about substructure on scales smaller than the size of 
the grid cells is ign ored. As a cure, adaptive grids 
have been invoked (e.g. Razoumov fc Sommer-Larsenll2006l . 
lAlvarez. Bromm. fc Shapiro! 120061). In SPH, grid-less ray- 
tracin g has been introduced by iKessel-Devnet fc Burkert] 
(2000) to investigate the effects of ionizing UV radiation by 
massive stars on the surrou nding interstellar medium. The 
Stromgren volume method (|Dale. Ercolano."fc Clarkell2007l , 
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compare also I Johnson. Greif. fc Brommll2007j) and the ion- 
ization front tracking technique employed by lYoshida et all 
(2007) are related approaches. Another grid-less radiation 
tracing scheme that has been app l ied to SPH simulations 
is presented in Ritzerveld fc Ickei (120061 ) . Inspired by the 



iKessel-Devnet fc BurkertJ (|200Ch approach. ISusal (|2006l ) de- 
scribes a radiation-hydrodynamics scheme for the transport 
of ionizing radiation in cosmological simulations of structure 
formation. 

One major drawback most ray-tracing schemes share 
with the Monte-Carlo technique is that the computational 
effort to solve the radiative transfer equation scales lin- 
early with the number of sources in the simulation box (but 
see Section T4.2.3I later on). For this reason, SPH radiation- 
hydrodynamical simulations employing the ray-tracing ap- 
proach have mostly been restricted to the study of problems 
involving only a few sources. 

In the following sections we will describe traphic, a 
novel method to solve the radiative transfer equation in 
SPH. traphic employs a photon tracing technique that 
works directly on the discrete set of irregularly distributed 
SPH particles. It is designed to overcome the challenges set 
by our goal of carrying out large radiation-hydrodynamics 
SPH simulations exhibiting a wide dynamic range in length 
scales and containing a large number of light sources. 



TRAPHIC 
CONES 



TRANSPORT OF PHOTONS IN 



In this section we give a detailed description of traphic, our 
method to solve the radiative transfer equation in SPH. We 
start the presentation of our radiative transfer scheme by 
sketching the main idea, in order to give an overview and to 
introduce the reader to the underlying concepts, which will 
be illustrated and explained in more detail in the subsequent 
sections, traphic can be applied to solve the radiative trans- 
fer equation in both two-dimensional and three-dimensional 
SPH simulations. When illustrating concepts with schematic 
figures we will, however, restrict ourselves to two dimensions 
for the sake of clarity. 

Although we ultimately aim at performing simulations 
in which the radiation transport is fully coupled to the hy- 
drodynamics, here we assume that the SPH density field is 
static and that the radiation does not exert any thermo- 
dynamic or hydrodynamic feedback on the matter. We will 
report on the coupling of radiation and hydrodynamics in 
future work. 



4.1 Overview 

We solve the radiative transfer equation (Eq. [7} in finite 
steps (t r ,t r + At r ), where t r is the current simulation time, 
by tracing photon packets radially away from their location 
of emission until a certain stopping criterion is fulfilled. We 
do not introduce a grid on which the photons are propa- 
gated. Instead, we propagate the photons directly on the 
discrete set of particles in the simulation. For the transport 
of photons we employ the same particle-to-neighbour com- 
munication scheme that is already used in the SPH simula- 
tion to solve the equations of hydrodynamics. That is, we 
accomplish the global transport of photons by propagating 



them only locally, between a particle and its Nngb neigh- 
bours. We allow Nngb to be different from N ng b, the number 
of neighbours used during the SPH calculations. In SPH, 
N ng b determines the adaptive spatial resolution of the hy- 
drodynamical simulation. Similarly, in our scheme, N ng b sets 
the adaptive spatial resolution of the radiation transport. It 
is the first of two numerical parameters in our method (see 
also the left-hand panel of Fig. [T} . 

Working directly on the set of SPH particles, our scheme 
exploits the full spatial resolution of the hydrodynamical 
simulation. A further advantage of our approach is that the 
radiation transport is automatically parallel on distributed 
memory, as long as the SPH simulation itself is parallel 
on distributed memory. This is in contrast to radiation- 
hydrodynamical schemes that employ individual communi- 
cation schemes for the transport of radiation and the SPH 
computations. Our radiative transfer scheme can hence be 
coupled to the SPH simulation without introducing any ad- 
ditional computational structures. 

The main challenge we face when performing radiation 
transport by propagating photons only from a particle to its 
neighbours is achieving a transport that is directed: Pho- 
tons travel along straight rays, while the spatial distribu- 
tion of the SPH particles is generally highly irregular. In 
the following we give an overview of the concepts employed 
in our scheme and describe how we overcome this and other 
challenges to accomplish the transport of radiation in SPH 
simulations. 

We distinguish two types of particles present in SPH 
simulations that are relevant for the transport of photons: 
Star particles and SPH, or gas, particles. Only gas particles 
can interact with photons, via absorptions and scatterings. 
We assume, however, that both star and gas particles can be 
sources of photons and will therefore refer to them as source 
particles, or simply sources. 

The transport of radiation at simulation time t r over a 
single radiative transfer time step At r starts with the emis- 
sion of photons by the source particles. Subsequently, these 
photons are propagated downstream, radially away from the 
source, simultaneously for all sources. In our particle-to- 
neighbour transport scheme photons are propagated from 
gas particle to gas particle. On their way, photons may in- 
teract with the matter field discretized by the gas particles. 
Each gas particle can simultaneously receive photons, possi- 
bly emitted at different times, from multiple (in principle all) 
sources in the simulation. Their further propagation would 
require a loop over all propagation directions. We explicitly 
avoid the resulting linear scaling of the computation effort 
with the total number of sources by introducing a source 
merging procedure. 

The distribution of photons amongst the neighbours of 
the source particles must respect the angular dependence 
of the source emissivity. We achieve this by introducing for 
each source a set of randomly oriented, tessellating emission 
cone f](SectionE23), as illustrated in the middle panel of 
Fig. [T] The fraction of the total number of photons emitted 
into each cone is proportional to its solid angle. The emission 



1 We use the word cone in a general sense. It does not necessarily 
describe, but includes in its meaning, a regular cone, which is 
defined as a pyramid with a circular cross section. 
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direction (which we will equivalently refer to as the propaga- 
tion direction) associated with each emitted photon packet 
is determined by the central axis of the emission cone. The 
number of cones used in the emission cone tessellation, N c , 
is the second numerical parameter in our method. It sets 
the angular resolution. The random orientation given to the 
cones ensures that photons are not restricted to propagate 
along a fixed number of directions and prevents artefacts 
introduced by the shape of the emission cones. 

Some of the emission cones may not contain any neigh- 
bours, as is the case for the bottom-right cone of the tes- 
sellation shown in the middle panel of Fig. [T] To neverthe- 
less emit photons into the corresponding directions, we cre- 
ate additional neighbours to which the photons can then 
be transferred (Fig. [1] right-hand panel) . We refer to these 
neighbours as virtual particles ( ViPs) , since they do not take 
part in the SPH simulation. The properties of the ViPs 
are obtained from those of the neighbouring SPH parti- 
cles using SPH interpolation. ViPs compensate for the lack 
of neighbours in certain solid angles around the emitting 
source. Hence, by employing ViPs we introduce the freedom 
of choosing emission directions independently of the geom- 
etry of the SPH particle distribution. As we will argue in 
Section 14.2.11 and investigate in more detail in Appendix 
|A1 this freedom is a necessary requirement for achieving a 
desired angular dependence of the transport process, e.g. 
the isotropic emission of photons by star particles, in any 
particle-to-neighbour transport scheme. 

The photon packets distributed amongst the neighbours 
of the sources are further transferred downstream until they 
experience an interaction. In contrast to the emission of 
photons by source particles, gas particles distribute photons 
only amongst the subset of neighbours located in the down- 
stream direction (see Fig. [2}. We distinguish this directed 
particle-to-neighbour transport from the emission process 
by referring to it as transmission (Section I4.2.2[l . In more 
detail, each gas particle transmits photons downstream by 
distributing photon packets amongst the subset of neigh- 
bours located in certain regular cones only. These so-called 
transmission cones are centred on the emission (propaga- 
tion) direction associated with each photon packet and have 
an opening angle determined by the angular resolution, i.e. 
by N c . They confine the photon packets to the cone into 
which they were emitted by the corresponding source parti- 
cle. As a result, the photon packets are propagated radially 
away from the sources, in a manner that is independent of 
the geometry of the SPH particle distribution. Like emission 
cones, transmission cones might not contain any neighbours. 
Again we employ ViPs to accomplish the photon transport 
into the corresponding directions. 

By keeping the opening angle of the cones fixed, the 
solid angle subtended by a transmission cone as viewed 
from the original source decreases with the distance from 
the source. As a result, the photon transport becomes adap- 
tive in angle, as in the ra y-splitting technique em ployed in 
ray-tracing schemes (e.g. lAbel fc Wandelj |2002| ). We will 
demonstrate later on, when testing a numerical implementa- 
tion of our scheme fSection !5.3.2|l . that the use of implicitly 
adaptive transmission cones confining the photon propaga- 
tion yields sharp shadows behind opaque obstacles. 

Gas particles receive photon packets from other par- 
ticles through the processes of emission and transmission 



described above. A given gas particle can simultaneously re- 
ceive multiple photon packets, which may each be associated 
with a different emission direction. In fact, our assumption 
that all star and gas particles in the simulation box can be 
source particles would imply that the number of directions 
from which photon packets can be simultaneously received 
by a given gas particle would scale linearly with the total 
number of particles in the simulation. The computational ef- 
fort to accomplish the subsequent photon transmission along 
the associated directions would then also scale linearly with 
the total number of particles. Since this consideration ap- 
plies to the transmission of photons by any SPH particle in 
the simulation, the total computational effort to solve the 
radiative transfer equation would scale no less than quadrat- 
ically with the total number of SPH particles. 

To avoid this computationally expensive scaling, we in- 
troduce a source merging procedure (cp. Fig. [4] which will 
be explained in more detail in Section|333J • We make use of 
the fact that source particles that are seen as close (in angle) 
to each other from a given gas particle can be approximated 
by, or merged into, a single point source. We average the 
emission directions associated with the photon packets re- 
ceived from sources in solid angle bins defined by a so-called 
reception cone tessellation which, except for a random ro- 
tation, is identical to the cone tessellation employed for the 
emission process. Consequently, photon packets need to be 
transmitted into at most N c directions. 

We distinguish two types of interactions with the mat- 
ter field sampled by gas particles and ViPs that photons 
may experience. Absorption interactions change the num- 
ber of photons contained in each packet. Scattering interac- 
tions change the propagation direction (and possibly the fre- 
quency) of the photons in a packet. All interactions are taken 
into account in an explicitly photon-conserving manner. We 
give a detailed description of interactions in Section T4. 31 

4.2 Transport of photons 

We now expand on the description of our radiative trans- 
fer scheme given above. We comment in detail on the main 
concepts underlying traphic, i.e. the emission of photons by 
source particles, the transmission of photons by gas particles 
and the source merging procedure. The description of how 
these concepts are employed to advance the solution of the 
radiative transfer equation in time is deferred to Section T4. 41 

4-2.1 Photon emission by source particles 

In this section we describe the emission of photons by source 
particles. Star particles emit photons according to their in- 
trinsic luminosity, which can for instance be obtained from 
population synthesis models. An example of the emission of 
photons by gas particles is the emission of recombination 
radiation by a recombining ion. 

Let us consider the emission of photons by source parti- 
cle i, located at r^. We denote the number of photons emit- 
ted per unit time per unit frequency v per unit solid angle 
Q around the unit direction vector n by A/L.,i(n, r, t), or sim- 
ply Afv,n,i- With this notation, the number of photons AC,i 
emitted per unit time at frequency v is J\f u ,i = f dQ Av,n,i, 
and the total number of photons emitted per unit time is 
A/" 7 ,i = / dv f dQ. Af v , n ,i- 
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Figure 1. Emission of photon packets by source particles. Left-hand panel: A source particle (black star) and its neighbourhood (big grey 
disc). In this example the radius h of the neighbourhood is defined such that there are N ngb = 4 neighbouring gas particles (white discs). 
Note that SPH simulations typically use N ngb 48. Middle panel: The randomly oriented emission cone tessellation of the source particle 
(solid lines) is shown for the case N c = 4. The dashed arrows indicate the central axes of the cones. Note that the bottom-right cone 
does not contain a neighbouring gas particle. Right-hand panel: The source particle has transferred its photon packets to its neighbours 
(black discs). The emission direction (small solid arrows) associated with each packet points in the direction of the axis of the cone in 
which the neighbour resides. For the empty cone, a virtual particle is created (red hatched disc), placed randomly along the central axis 
of the cone. 



Source particle i emits J^ +Atr _A/ 7 ^ photons over the 
radiative transfer time step At r - In agreement with our 
particle-to-neighbour based transport approach, the emis- 
sion process is modelled by distributing these photons 
amongst the N ngb nearest gas particles, residing in a sphere 
of radius hi centred on r^. This sphere is schematically de- 
picted in the left-hand panel of Fig. [T] 

Source particle i distributes its photons amongst its 
neighbours with the help of a set of space-filling emission 
cones, defined as follows (middle panel of Fig. [TJ . We tes- 
sellate the simulation box into N c cones of (generally not 
identical) solid angles Of, k — 1,2, N c , with the apexes 
located at the position of particle i. The number of neigh- 
bours N* gbi in each cone k is determined, taking values in 
the range < N* gb>i < N ngb . 

For what follows we consider the emission of photons 
for each frequency v separately. To each neighbour j in cone 
k a photon packet of characteristic frequency v is emitted. 
The packet contains a fraction > (j = l...N^ gbi ) of 
the total number of photons N v .i to be emitted per unit time 
at frequency v, with 

w i = —f = x — . (8) 

J dQ Mv,n,i y"ngb,i w i 



where w\ are weights attached to neighbour j in cone k. The 
first factor on the right-hand side of Eq.[8]determines which 
fraction of the total number of photons is emitted into cone 
k, whereas the second factor controls the fraction of photons 
that is transferred to neighbour j, that is, it controls the 
distribution of photons amongst the particles within cone k. 
Here we set w\ = 1, i.e. equal weights for all neighbours in 
a given cone. The number of cones used in the tessellation, 
the parameter N c , determines the angular resolution of the 
radiation transport. A cone tessellation may consist of cones 



with different solid ang We therefore define the angular 
resolution to be the size of the average solid angle, (f2) = 
4n/N c . 

For each emission cone k the central axis is defined, 
characterised by the unit vector pointing away from the 
source (Fig. [T] middle panel). We refer to this vector as 
the emission direction associated to photon packets emitted 
into cone k. When a photon packet is emitted to a neigh- 
bouring gas particle, the emission direction is transferred in 
addition to the number of photons it contains (Fig. [1] right- 
hand panel). Since the orientation of the emission cone tes- 
sellation is randomised by applying a random rotatioijfl at 
each emission, there is no a priori limit on the values the 
emission directions can take. Note that while we transfer 
the emission direction, we do not transfer the position of 
the source. Photon packets are traced further downstream 
based on their emission direction only, as we will explain 
in Section \4. 2. 2 1 Each photon packet has also an associated 
clock t*. At emission the clock time is set to the time of 
the simulation, t* — t r . In Section [4.41 we will see that the 
clock can be employed to propagate the photon packets at 
the speed of light. 

Some cones may not contain any neighbouring gas par- 
ticles. This is for instance the case for the bottom right cone 
in the middle panel of Fig. [1] For a fixed number of neigh- 
bours these empty cones will occur more frequently if the 
angular resolution is higher (i.e. if the solid angles of the 
cones are smaller). On the other hand, for a fixed angular- 
resolution N c , empty cones will occur more frequently if the 
number of neighbours N ngb is smaller. Spatial clustering of 
the neighbours also increases the probability of the occur- 
rence of empty emission cones. In the absence of a neigh- 

2 We note that in the numerical implementation of our radiative 
transfer scheme the tessellation is made of cones having equal 
solid angles, see Appendix [Bl] 

3 For a numerical implementation of the random rotations we 
refer to Appendix IB2I 
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bouring gas particle photons cannot be transported along 
the emission direction of the empty cone. We therefore cre- 
ate a new neighbour, a so-called virtual particle (ViP), to 
which the photons are transferred. The ViP is placed along 
the cone axis at a random distance < hi from the source par- 
ticle, such that the volume of the cone is uniformly sampled 
(see the right-hand panel of Fig. [1]). The properties of ViPs 
(e.g. density) are determined from their neighbouring gas 
particles using SPH interpolation. ViPs are only employed 
for the transport of photons. We stress that ViPs are not 
used by the SPH simulation. 

We emphasize that the introduction of cones is essen- 
tial for accomplishing the emission process in our particle- 
to-neighbour transport scheme. To see this, consider the al- 
ternative of distributing the photons directly amongst the 
neighbours of gas particle i. This amounts to setting Qj — 4-k 
and 0,^ = for j > 1 in Eq. [8] Without any reference sys- 
tem, the emission directions associated to the photon pack- 
ets would then be given by the unit vectors pointing from 
source i to neighbour j. In Appendix lAl we show that in 
this case the net emission direction in general would cor- 
relate with the direction towards the centre of mass of the 
neighbours. As a result, the emission process would depend 
on the clustering of the neighbours in space as set by the 
geometry of the SPH simulation. The use of emission cones 
combined with ViPs gives us the freedom to choose emis- 
sion directions independently of the spatial clustering of the 
neighbours. This allows us to model any angular dependence 
of the source emissivity, within the bounds of the chosen an- 
gular resolution. 

In summary, source particles transfer photon packets 
to their neighbouring gas particles using a set of randomly 
oriented, tessellating cones. Each cone defines a direction 
of transport, i.e. the emission direction, associated with the 
packets. The random orientation applied to the cone tes- 
sellation ensures that photon packets are not transferred 
along a fixed set of directions only. Virtual particles (ViPs) 
are placed in emission cones not containing any neighbours, 
to which the photon packets are then transferred to. The 
combination of emission cones and ViPs makes the radia- 
tion transport independent of the spatial distribution of the 
neighbours. 

4-2.2 Photon transmission by gas particles and ViPs 

In the last section we described how a source particle dis- 
tributes photon packets amongst the gas particles in its 
neighbourhood. For each packet we defined an emission 
direction, independently of the spatial distribution of the 
neighbouring gas particles by employing a randomly ori- 
ented set of emission cones. In this section we describe how 
the packets are propagated through the simulation box along 
their emission directions, by employing directed particle-to- 
neighbour transport, a process which we refer to as trans- 
mission. 

Consider a gas particle i, which receives a photon 
packet. Recall that the neighbours of particle i are the N ng b 
nearest gas particles located in the sphere with radius hi, 
centred on particle i (Fig. [2] left-hand panel). Analogous 
to the case of emission, particle i re-distributes the photons 
contained in the packet amongst its neighbours. In contrast 
to emission, photons are only distributed amongst the sub- 
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Figure 3. The apex angle ui of the transmission cones as a func- 
tion of the angular resolution N c (see Eq. |9)l . 

set of neighbours N t , n gb ^ N ng b located downstream. These 
are the neighbours residing in a regula^ cone centred on the 
direction of propagation of the packet (see the middle panel 
of Fig. [2j. We refer to this cone as transmission cone. 

The apex of the transmission cone is located at the po- 
sition of gas particle i. The solid angle of the cone is set by 
the angular resolution, Q t = 4n/N c = (Q), and determines 
the apex angle u through the standard relation 



i . Ot\ 180 

u! = 2 arccos 1 x . 

2ty 7T 



(9) 



We show this relation in Fig. [3] By definition, a neighbour 
j with position rj is interior to the transmission cone with 
apex at the position n of the transmitting particle i if the in- 
ner angle between the transmission cone axis and the vector 
rj — r» is less than ui/2. 

The received photon packet is split into Nt, ng b packets, 
each of which is transferred to one of the downstream neigh- 
bours j of particle i. The number of photons contained in 
each packet is set by the weights «r? , such that each neigh- 
bour j receives a photon fraction wl / ~^2 k wf of the parent 
photon packet, where the sum is over all downstream neigh- 
bours k = 1, Nt,ngb- Here we assume wj = 1, giving equal 
weight to each neighbour. If there are no downstream neigh- 
bours, i.e. if the transmission cone is empty, we simply create 
a neighbour as in the case of empty emission cones described 
in the previous section. That is, we place a virtual particle 
(ViP) along the emission direction of the photon packet, 
a random (in volume) radial distance < hi away from the 
transmitting gas particle. The properties of the ViP are de- 
termined by SPH interpolation from its neighbours, and the 
photon packet is transferred. 

Each photon packet inherits the emission direction from 
its parent packet (Fig. [2] right-hand panel). After all pack- 
ets have been transferred to the downstream neighbours, the 



A regular cone is a pyramid with a circular cross-section. 
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Figure 2. Transmission of photon packets by gas and virtual particles. The particle positions shown are the same as in Fig. [T] as are 
the neighbourhood (dashed circle) and emission cone tessellation of the source particle (white star). The transmission of photons is 
illustrated for one of the neighbours that received radiation from the star in Fig. [T] Left-hand panel: The neighbourhood (big grey disc) 
of the transmitting gas particle (black disc) is defined. As in Fig. [T] N ng b = 4. The emission direction associated with the photon packet 
to be transmitted is shown as the short solid arrow. Middle panel: The transmission cone is shown, centred on the emission direction of 
the photon packet that is to be transmitted. In this cone one downstream neighbour is found. Right-hand panel: The photon packet is 
transmitted to the downstream neighbour, turning it into a transmitting particle. The cycle repeats with defining the neighbourhood for 
this particle. As a result, the photon packet is propagated downstream, radially away from the source particle. 



transmission is finished. Subsequently, the transmission pro- 
cedure can be applied again. The transmission process thus 
generally splits each photon packet into multiple packets, 
propagating with the same emission direction. With each 
subsequent transmission individual photon packets are con- 
fined to a smaller fraction of the solid angle traced out by 
the emission cone they were emitted into (see Fig. middle 
panel) . This is similar to the technique of ray splitting em- 
ployed in ray-tracing codes. Hence, in traphic the photon 
transport takes place adaptively in angle. 

The transmission procedure specified above for gas par- 
ticles is also applied for ViPs. Again, the transmission cone 
of the ViP can either contain gas neighbours or be empty. 
If it is empty, the ViP creates another ViP, as in the case of 
transmission by gas particles. After the ViP has performed 
the transmission, it is deleted. ViPs are therefore temporary 
constructs. For N c N ng b the total number of ViPs in the 
simulation is proportional to N c , whereas for N n gb S> N c 
the total number of ViPs approaches zero. 

The main purpose of the transmission cones is to con- 
fine the downstream propagation of photons to the solid 
angles into which they were emitted by the source particles, 
which is the main challenge for schemes using particle-to- 
neighbour transport. Just as emission cones were introduced 
to make the emission of photons by source particles inde- 
pendent of the geometry of the SPH particle distribution, 
transmission cones are introduced to further propagate the 
photons downstream, independently of the geometry of the 
SPH particle distribution. As a consequence, shadows will 
be thrown behind opaque obstacles, and their sharpness will 
be in agreement with the chosen angular resolution. In fact, 
as we will see in Section lo.3.21 the shadows are much sharper 
than implied by the formal angular resolution. 

The cones making up the emission tessellation will gen- 
erally be irregular - in three dimensions there exists no tes- 
sellation of space with regular cones. The use of regular cones 
for the downstream propagation of photon packets emitted 
into a certain irregular emission cone is therefore an approxi- 



mation. In principle, each photon packet could be transmit- 
ted into a cone having the shape of the emission cone it 
originates from. However, tracing photons within irregular 
cones is computationally more demanding. Furthermore, in 
the next section we will see the necessity for combining mul- 
tiple photon packets propagating along different directions 
into a single photon packet that propagates along a corre- 
spondingly averaged direction. Since there is no unique way 
of defining the shape of a transmission cone associated with 
this average direction, we use a regular shape. 



4-2.3 Photon merging by gas particles 

In the previous section we described the transmission of a 
single photon packet arriving at a gas particle. In principle, 
each ga^f] particle can simultaneously receive multiple pho- 
ton packets, possibly emitted by different sources at different 
times, propagating along different emission directions, a sit- 
uation which is depicted schematically in the left-hand panel 
of Fig. 3] Accomplishing the transmission would then require 
executing a loop over all received packets for each transmit- 
ting gas particle. Clearly, for each gas particle this would 
imply a linear scaling of the computation effort with the 
number of packets that need to be transmitted, or in other 
words, with the total number of source particles present in 
the simulation. 

In cosmological simulations of structure formation, for 
instance, where a significant number of star particles can 
usually be found already at high redshift, the simultaneous 
reception of multiple photon packets will be the rule rather 
than the exception. Even without taking into account the 
diffuse radiation of the intergalactic gas resulting from radia- 
tive de-excitations, recombinations and scatterings, a linear 
scaling with the number of sources would quickly turn the 



5 We do not consider ViPs here, since by construction each ViP 
can receive only a single photon packet. 
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Figure 4. Merging of photon packets by gas particles. Left-hand panel: A gas particle (white disc) is receiving photon packets simulta- 
neously from two transmitting gas particles (black discs). The emission direction associated with each photon packet is indicated by an 
arrow, whose length is given by the number of photons contained in the packet. The big grey discs indicate the neighbourhoods of the 
transmitting gas particles, the solid lines are the transmission cones (we assume N c = 4). For clarity, we do not show particles and photon 
packets that are not directly involved in the merging event. Middle panel: A randomly oriented reception cone tessellation is defined for 
the receiving gas particle. We omitted the transmission cones and the neighbourhoods of the transmitting particles. Right-hand panel: 
The photon packets have been transferred to the receiving gas particle, turning it into a transmitter. Their emission vectors (grey arrows), 
after translation to the position of the receiver, are both found to fall into the upper-right reception cone. Note that the particles itself fall 
into different reception cones. Their positions are irrelevant for the merging, since transmitted photon packets rather than transmitting 
particles are merged. The emission vectors are merged through a weighted average as described in Section 14.2.31 resulting in a single 
photon packet (black arrow). 



radiation transport into a computational task that would be 
too expensive to be carried out even with the most advanced 
supercomputer available today. 

A few radiation tracing schemes have attempted to 
tackle this scaling problem by replacing source s that are 



close to each other with a single point source (e.g.lCenll2002l: 
iRazoumov fc Sommer-Larseni 1200(1 iGnedin fc Abell l200lf ). 
Here we introduce a photon packet merging procedure that 
strictly limits the number of packets that need to be trans- 
mitted per particle to at most N c , without restricting the 
number of directions along which each individual packet can 
propagate. 

For each gas particle i receiving a photon packet we de- 
fine a so-called reception cone tessellation, as shown in the 
middle panel of Fig. [4] A reception cone tessellation is a 
tessellating set of N c cones attached to a gas particle. It is 
identical to the set of N c cones employed for the emission 
of photons by source particles, but generally has a different 
orientation. As their name indicates, and in contrast to the 
emission cones, reception cones are used to collect photons. 
Whenever a photon packet p is transferred to gas particle 
i, the packet is binned into one of the reception cones by 
examining into which reception cone its emission direction 
n p falls, after translation to the location of the gas particle. 
Photon packets whose emission directions fall in one and 
the same cone are merged, leaving only a single packet for 
transmission (Fig. [4] right-hand panel). The reception cone 
tessellation is given a random orientation to prevent arte- 
facts. 

A merged photon packet created in cone k is assigned 
the new, merged emission direction n m ,fc. This new emis- 
sion direction is defined as the weighted average n mj fc = 
^2 Wprip/ ^2 Wp, where the sum is over all photon packets re- 
ceived in emission cone k (k = 1, 2, N c ), and the weights 
w p are the number of photons per packet p that are to be 



transmitted into direction n p to the downstream neighbours 
of particle i (Fig. 01 right-hand panel). Similarly, the merged 
photon packet is associated a merged clock t* m k by averag- 
ing over the individual clock times tp~ in reception cone k, 

i-e- = E w p^/E«V 

As a result of the photon packet (resp. source) merging 
at most N c packets have to be transmitted per gas particle, 
fully consistent with the angular resolution of our trans- 
port scheme. The computation effort required to perform 
radiative transfer simulations with traphic can therefore 
be readily controlled: In the most demanding situation, i.e. 
in the optically thin limit, when the whole box is filled with 
radiation, it merely scales with the product of spatial and 
angular resolution, N x N ng i> x N c , where N is the total num- 
ber of particles (SPH + stars). The merged photon packets 
are not associated with any existing source particle in the 
simulation. They should be thought of as emerging from an 
imaginary source, whose properties are implicitly defined by 
our merging prescription. The merged photon packets are 
traced further downstream, radially away from this imagi- 
nary source, into the merged emission direction, exactly as 
was described for non-merged emission directions in the pre- 
vious section. 



4.3 Photon interactions with gas particles 

Photons are propagated radially away from each source until 
they interact with the matter represented by the SPH par- 
ticles. Here we distinguish between absorptions and scatter- 
ing interactions and describe how they are accounted for in 
traphic. We remind the reader that in this work we do not 
discuss the effect of interactions on the thermodynamical 
and hydrodynamical evolution of the SPH particles. 
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4-3.1 Absorption 

Absorptions remove photons from the beam emitted by 
the source. This process is described in terms of the 
mass absorption coefficient k„. From the formal solution of 
the radiative transfer equation it can easily be seen (e.g. 
iMihalas fc Weibel Mihalasl GUI) that the fraction of pho- 
tons of frequency v that is absorbed while travelling along 
the straight line connecting the positions ri and r2 is given 
by 1 — exp(— r), where 

r = f 2 dr K v {r)p{r), (10) 

is the optical depth over the distance d = \r 1 — r2 1 . 

In traphic, absorptions are accounted for by remov- 
ing the photon fraction 1 — exp( — ry) from photon packets 
propagating between particle i and its neighbouring parti- 
cles j. For the calculation of the optical depth Tij between 
these two particles we need to numerically evaluate the in- 
tegral Eq. 1101 This is computationally expensive, since it 
involves the calculation of the density p (and the absorption 
coefficient k„) at a large number of points along the photon 
path using the SPH formulation . Similar to the approach of 
iKessel-Devnet fc Burkertl (|2000T ) . we therefore approximate 
the density field near a gas or virtual particle i by the SPH 
density pi evaluated at its position. 

Since photon packets are propagated between particle i, 
located at r^, and particle j, located at Tj, along their emis- 
sion direction, i.e. in the direction of the unit vector n, the 
distance d they cover is generally not equal to the distance 
between the particles. Instead, we set d = |n- (r^ — r 3 -)[, i.e., 
we employ the projection of the particle distance along the 
emission direction. This is a valid approximation far away 
from the source from which the photon packets originate, 
but it generally fails close to the source. Therefore, we only 
employ the projected distance when propagating photons 
between a transmitting particle and its neighbours. On the 
other hand, when propagating photons from a source parti- 
cle to its neighbours, we use the particle distance, i.e. we set 
d = \ij — r»|, corresponding to radially outward propagation. 

The distance d is used to obtain the optical depth Ty 
between particle i and particle j, Tij — n v pjd. This approx- 
imation neglects the existence of any substructure between 
particle i and its neighbouring particle j. However, this does 
not result in a loss of information if the radius hi of the 
neighbourhood of particle i is chosen such that hi < hi 
(resp. N ngb < N ngb ). 

The distance d is furthermore used to advance the clock 
t* of each photon packet according to the rule 

t*^t* + d/c. (11) 

By synchronising the clock t* with the simulation time t r , 
the clock can be employed to advance photon packets at 
the speed of light, as we discuss in Sec. 14.41 We note in 
passing that the clock t* can also be used to implement the 
cosmological redshifting of photons. 

The number of absorbed and transmitted photons is 
now calculated as follows. From the photon packet emit- 
ted or transmitted by particle i to particle j, a fraction 
1 — exp( — Tij) is absorbed by particle j. The photon packet is 
subsequently transmitted by particle j containing a fraction 



exp(— Tij) of the original number of photons. By construc- 
tion, this procedure explicitly conserves photons, since 

1 - e~ Ti ' + e~ Til = 1. (12) 

Photons absorbed by ViPs are distributed amongst all 
neighbours of the ViP using (photon-conserving) SPH inter- 
polation. This is necessary, because ViPs are a temporary 
construct employed for the transport of photons; permanent 
information is only stored at the particles in the SPH simu- 
lation. For further reference, we denote the total number of 
photons impinging on and the total number of photons ab- 
sorbed by particle i over a single time step At r with AA/in,i 
and AAf a bs,i, resp. 

4.3.2 Scattering 

In contrast to absorptions, scatterings change the direction 
and possibly the frequency of the interacting photons. A 
scattering event can be thought of as an absorption event 
followed by an immediate re-emission. Scatterings can hence 
be described by Eq. [7J after adapting the absorption coeffi- 
cient k„ and the emissivity e„ . The re-emission re-distributes 
the photons in angle (and possibly frequency). For this rea- 
son the scattered photons are sometimes referred to as dif- 
fuse. Adapting the emissivity for scatterings requires eval- 
uating an integral over the intensity. Therefore, scatterings 
turn the radiative transfer equation (Eq. [7J into an integro- 
differential equation. 

In traphic, scatterings, that is the transport of diffuse 
photons, are modelled by a combination of absorption and 
re-emission events. Photons to be scattered travelling be- 
tween two particles are first removed from the photon packet 
as described in the last section. For a true absorption event, 
the photon energy would be lost to the thermal bath pro- 
vided by the matter. In contrast, in case of scattering, a gas 
particle that absorbed photons re-emits them, closely fol- 
lowing the description for source particles in Section [4.2.11 

In particular, we employ randomly oriented emission 
cones in combination with Eq. [5] to re-emit the absorbed 
photons. Since modelling the scattering process means fol- 
lowing photons emitted earlier by a source, the clock of the 
photon packets that are to be re-emitted is not set to the 
simulation time t r , as was the case for the intrinsic emission 
of photons by source particles. Instead, it is determined by 
the clock of the photon packets that are scattered and the 
details of the scattering interactions (e.g. there could be a 
time delay between absorption and re-emission, or photons 
could have travelled a distance greater than the (projected) 
particle distance d when inferring the properties of the scat- 
tering process from a sub-resolution model) . If the scattering 
details are taken into account by correspondingly adjusting 
d, Eq. [TT]can be employed to find the new clock time. Note 
that it is crucial for the modelling of scatterings to employ 
a radiative transfer scheme which does not scale with the 
number of sources, since every gas particle will soon re-emit 
photons. 

4.4 Solving the radiative transfer equation 

traphic solves the radiative transfer equation (Eq. [Jj in 
discrete time steps At,., the size of which depends on the 
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details of the problem under study and has to be decided on 
a case by case basis. Hence, we defer its discussion to the 
second part of this work, where we present an implementa- 
tion of our method to solve specific problems. In this section 
we briefly review the concepts described in the previous sec- 
tions used to obtain the solution to the radiative transfer 
equation over a single time step before discussing how this 
solution is advanced in time. 

Starting at simulation time t r , the size of the radiative 
transfer time step At r is determined. Source particles emit 
photon packets to their N„ g b neighbouring gas particles em- 
ploying N c cones. The total number of photons contained 
in the packets is determined by the number of photons N v ,i 
to be emitted per unit time in frequency bin v and by the 
size of the time step At r . For each neighbouring gas parti- 
cle the optical depth to the source particle is evaluated and 
the number of photons interacting with the matter on their 
way from the source particle is inferred. Some photons may 
be absorbed, some may be scattered; the remaining photons 
are left unimpeded for transmission. 

For the subsequent transmission of photons by the gas 
particles, photons associated to sources that are seen as 
close (in angle) to each other are merged, limiting the num- 
ber of directions into which photons have to be transmitted 
to at most N c . Next, for each transmitting and scattering 
gas particle the N ng b neighbours are found. Photons to be 
transmitted are distributed amongst the subset N t ,„ g b of 
neighbours that are located downstream, as determined by 
the corresponding (merged) emission directions. Photons to 
be scattered are distributed amongst all N ng b neighbours 
following the prescription in Section 14.31 Virtual particles 
(ViPs), which are placed in cones that do not contain any 
neighbours, are deleted after having transmitted or scat- 
tered their photons. Photons that were absorbed by a ViP 
are distributed amongst the neighbours of the ViP using 
(photon-conserving) SPH interpolation. 

The cycle described in the previous paragraph repeats 
until a stopping criterion is satisfied. The form of the stop- 
ping criterion depends on whether one solves the time- 
independent or the time- dependent radiative transfer equa- 
tion. When solving the time-independent radiative transfer 
equation, i.e. when neglecting the first term on the left hand 
side of Eq.[7] one formally assumes c — > oo. Accordingly, the 
stopping criterion becomes independent of the speed of light. 
The cycle can then for example be repeated until all pho- 
tons have been absorbed or have left the simulation box. In 
contrast, when solving the time-dependent radiative trans- 
fer equation, the stopping criterion is directly tied to the 
speed of light c: the cycle is repeated until all photons have 
propagated a distance cAt r . 

To decide whether or not a photon packet has travelled 
over the distance cAt r , we employ its clock t*. For further 
illustration it is convenient to explicitly follow the clock of 
a photon packet as it is propagated through the simulation 
box. Upon emission, the clock of the photon packet is ini- 
tialized with the simulation time t r , as already mentioned in 
Sec. 14.2.11 After the emission, when the photon packet has 
been transferred from the source to one of the neighbouring 
SPH particles (or to a ViP), its clock is advanced according 
to Eq. [11] It is then checked whether the clock has reached 
or crossed the threshold value 



t* th = t r + At r . (13) 

If not, the photon packet is transmitted further downstream. 
As before, its clock is advanced according to Eq. 1111 The 
photon packet is repeatedly transmitted and its clock is ad- 
vanced until it has reached or crossed the threshold value 
defined in Eq. 1131 At this point, the value displayed by its 
clock can in full generality be expressed as t* = t* th + eo, 
with eo ^ 0. 

The propagation error, eo, is typically larger than zero 
because the set of particles on which the photon packets are 
propagated is discrete, with the particle-to-neighbour dis- 
tances generally not related to cAt r . Employing Eq. 1131 to 
stop photon packets therefore results in the photon packets 
typically being propagated too far. Due to the Lagrangian 
nature of the SPH simulation, eo is individual for each pho- 
ton packet. If d is the particle-to-neighbour distance corre- 
sponding to the emission or transmission event after which 
the photon packet was stopped, then eo < d/c. The propa- 
gation error eo is therefore strictly limited, eo < maxi hi/c, 
where the maximum is over all emitting and transmitting 
particles. Note that eo becomes smaller with higher spatial 
resolution. 

Once stopped, a photon packet is held (and hence its 
clock stands still) until the simulation time progresses to a 
value t r + At r > t* , at which point the packet is transmitted 
and its clock is advanced again. The clocks of photon packets 
are thus kept in synchronization with the current simulation 
time, which effectively matches their propagation speed to 
become the speed of ligh10. Photon packets may be held over 
several radiative transfer time steps, as eo may be (several 
times) larger than At r . While the photon packet is held, 
its propagation error with respect to the current simulation 
time steadily decreases with each passing radiative transfer 
time step. Since t* is fixed while t r is increasing in steps of 
At r , the current propagation error decreases according to 
e n = e n _i — Atr, where n indicates the number of radia- 
tive transfer time steps that have passed since the photon 
packet was stopped. Immediately before the photon packet 
is propagated again, its propagation error is strictly limited 
by the size of the time step, — At r ^ e n < 0. The photon 
packet is then propagated such as to again bring the clock 
into synchronization with the simulation time, as described 
above. 

When multiple photon packets are merged into a single 
packet, the clock of the merged packet is determined by the 
averaging procedure described in Section [4.2.31 As a result, 
photons may be propagated over distances that differ some- 
what from the case without merging. We have seen above 
that, due to the synchronization of the photon packet prop- 
agation with the current simulation time t r , clocks display 
only values in the interval t r ^ t* ^ t r + At r + eo. The differ- 
ence in the clocks of different photon packets and hence the 
photon propagation error introduced by the merging proce- 
dure can therefore be controlled to become arbitrarily small 
by increasing the temporal (i.e. decreasing At r ) and spatial 
(leading to smaller eo) resolution. 

6 We note that the clocks can also be employed to propagate pho- 
ton packets at speeds c different from the speed of light (cp. the 
reduc ed speed of light approximation suggested in lGnedin &: Abel 
2001), by simply making the replacement c — * c. 
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After the stopping criterion (time-independent or time- 
dependent) has been fulfilled for all photon packets, the 
properties of gas particles are updated according to the ra- 
diative interactions that occurred. We emphasize that the 
properties of the gas particles (e.g. the ionization state) are 
therefore only updated at the end of each time step. Within 
each time step the order in which photon packets from dif- 
ferent sources arrive at gas particles is therefore irrelevant. 
After updating the particle properties, the simulation time 
is advanced, t r —* t r + At r . Finally, the size of the next 
time step is determined and the radiation is transported as 
described above. 

In summary, we have presented a radiative transfer 
scheme for use in SPH simulations that works directly on 
the unstructured grid formed by the discrete set of irregu- 
larly distributed SPH particles, traphic thus employs the 
full spatial resolution of the SPH simulation. It achieves di- 
rected transport of radiation by adaptively tracing photon 
packets in cones, traphic can be used to solve both the 
time- independent and the time-dependent radiative transfer 
equation in an explicitly photon-conserving way. Our scheme 
is by construction parallel on distributed memory machines 
if the SPH simulation itself is parallel on distributed mem- 
ory machines. Furthermore, the computation time necessary 
to accomplish the radiation transport does not scale linearly 
with the number of sources. Instead, it merely scales with 
the product of spatial and angular resolution, making our 
scheme suitable for simulations containing a large number 
of sources as well as for taking into account a diffuse radia- 
tion component. 



4.5 Reduction of particle noise 

The radiative transfer equation describes the propagation of 
photons within a continuum medium. In our scheme, how- 
ever, photons are propagated on the discrete set of SPH par- 
ticles, localised at irregular positions dictated by the SPH 
simulation. Consequently, numerical noise may arise from 
the discreteness and irregularity of the spatial distribution of 
SPH particles and this could influence the numerical solution 
of the radiative transfer equation obtained with traphic. 

To see how this particle noise can be reduced, it is 
helpful to employ a formal analogy (e.g. iMonaghanl 120051 ) 
between estimating the density in SPH simulations and es- 
timating a probability density from sample points: We can 
consider the positions of the SPH particles as a randorrQ 
sample drawn from a probability density function propor- 
tional to the mass density. We can therefore Monte Carlo re- 
sampk-0 the density field without sacrificing its information 
content. Periodically assigning new positions to the SPH 
particles during the radiation transport by resampling the 
density field should lead to a better representation of the 
continuum physics by the discrete set of SPH particles and 
hence lead to a reduction of particle noise. 

7 We stress that SPH itself is no t a M onte Carlo method, as 
first noted in lGineold fc Monaghar] lll978T ) . We only appeal to the 
Monte Carlo picture for use with the radiative transfer, not for the 
hydrodynamic evolution of the particles in the SPH simulation. 

8 We explicitly note that the resampling does not affect the 
hydrodynamical simulation, for which the original positions are 
used. 



To employ the resampling, we think of particle i at po- 
sition n in the simulation box as being de-localised within 
its sphere of influence of radius hi centred on ri and as- 
sume that the probability of finding it in the volume d 3 r 
around a particular point r in that sphere is given by the 
value of the interpolation kernel at that point (cp. the scat- 
tering approach of Section [2]). For the radiative transfer on 
a static set of particles with positions we then period- 
ically Monte Carlo re-distribute the particles within their 
sphere of influence according to this probability. We note 
that while we change the positions of the SPH particles, we 
keep all their other properties (e.g. density) fixed in order 
to avoid introducing scatter in the physical variables. We 
describe the numerical implementation of the resampling in 
Appendix IB3I In Sections 15.3.11 and 15.3.21 we will demon- 
strate that the application of this recipe leads to a signifi- 
cant reduction of particle noise. We will, however, also see 
that for most applications the noise is small. We therefore 
do not employ the resampling technique in our simulations, 
unless explicitly stated. 

Since the resampling randomly offsets the SPH parti- 
cles from their positions provided by the hydrodynamical 
simulation, the apexes of the transmission cones attached 
to them are randomly offset, too. The resampling proce- 
dure may therefore lead to a slight diffusion of photons out 
of the emission cone they were emitted into, effectively de- 
creasing the angular resolution. Note that a similar diffusion 
of photons could occur when solving the radiative transfer 
equation coupled to the hydrodynamics (on which we will 
report in a future publication), because of the physical par- 
ticle motion. We will study this diffusion in our numerical 
implementation of traphic in Section f5. 3. 21 



5 APPLICATION - TRANSPORT OF 
IONIZING RADIATION 

In this section we apply the radiative transfer scheme pre- 
sented in Section [3] to the transport of ionizing radiation. 
Ionizing radiation is thought to play a key role in determin- 
ing the ionization state and shaping the spatial distribution 
of the baryonic matter in the universe on both small and 
large scales. Examples include the triggering and quenching 
of star formation through radiative feedbac k from nearby 
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and present-day universe (e.g. 
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), the thin shell instabil- 



and the growth and percolation of ionized regions generated 
by the first stars and quasars during the so-ca lled Epoch 
of Reionization (for recent simulations see e.g. Illiev et al.l 
2006al;lTrac fc Cenl2006l ; lKohler. Gnedin. fc Hamiltonll2007l ; 
07). 



Paschos et aill2007 

In the following we describe a numerical implementa- 
tion of traphic, our radiative transfer scheme, specified for 
the transport of mono-chromatic hydrogen-ionizi ng radia- 
tion, in the state-of-the-art SPH code gadget-2 (|Springell 
2005). We start in Section 15.11 with a short review of the 
photo-ionization process. In Section [5.21 we present the de- 
tails of our numerical implementation, establishing the con- 
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nection to the general description given in Section[4] Finally, 
in Section 15.31 we perform several radiative transfer simu- 
lations, comparing the performance of our numerical imple- 
mentation of traphic in three different test problems to 
both analytic solutions and numerical simulations carried 
out with other radiative transfer codes as reported in the 
literature. 



TionTree 
Tion T" Tr« 



(20) 



From Eq. \T§\ we see that the equilibrium ionized fraction 
is exponentially approached on the instantaneous ionization 
equilibrium time-scale r eq . We will employ this time-scale 
later on for the numerical integration of the rate equation. 



5.1 Photo-ionization rate equation 

Here we briefly recall the principles of the photo-ionization 
and recombination process occurring for a hydrogen-only gas 
of mass density p exposed to hydrogen-ionizing radiation. 
We will employ the equations derived in this section in the 
description of the numerical implementation of traphic. 

Hydrogen-ionizing photons can be absorbed by neutral 
hydrogen. We approximate the frequency-dependence of the 
mass abs orption coe fficient k„ for hydrogen-ionizing radia- 
tion fe.g. lOsterbrocl3 [l989h bv 



(To 



(14) 
(15) 



with rim = (l — x)p/ m H the neutral hydrogen number den- 
sity, vo the Lyman-limit frequency of energy h p vo = 13.6 eV, 
(Jo = 6.3 x 10 -18 cm 2 the absorption cross-section for pho- 
tons at the Lyman-limit, mu the mass of a hydrogen atom 
and <d(x) the Heaviside step function; the ionized fraction 
is x — nmi/TiH- The number of photo- ionizations per unit 
time per neutral hydrogen atom at a certain point in space 
is determined by the photo-ionization rate F, 



1 , 4nJ v (v) 

dv : (7„ 



(16) 
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where J„ = f dfi 7„/(47r) is the mean ionizing intensity. The 
rate of change of the neutral fraction 77 = 1 — \ a t this point 
is then 

± T] = a(T)n e x-rr,= - K -• (17) 

OX Tree Tion 

In the last equation, a(T)n e is the number of recombinations 
occurring per unit time per ionized hydrogen atom, r rec is 
the recombination time-scale and Ti 0n is the photo-ionization 
time-scale |f] 

With the definition \ = Tree/ (Tion+T r e C ) we can rewrite 
Ed. IT?lto read 

dX = X-X 

dt n on X ' 



(18) 



Setting dx/dt = yields the the equilibrium ionized fraction 
Xeq = Trec,eq/(Tion+T r ee,eq)- Over time-scales that are short 
compared with r rec /|dr rec /dt| and rie/\dn e /dt\, Eq. HSI con- 
stitutes a first order linear homogeneous differential equa- 
tion in x — X with constant coefficients, whose solution reads 

X(t)-Xe„ = (x(*o) - Xe,)e"^ (19) 

9 Note that in Eq. 1171 collisional ionizations can easily be taken 
into account by replacing T with (r + C(T)n e ), where C(T)n e 
describes the number of collisional ionizations per unit time per 
neutral hydrogen atom. In this work, however, we assume that 
collisional ionizations are unimportant, setting C = throughout. 



5.2 Numerical implementation 

We have adapted traphic for the transport of hydrogen- 
ionizing radiation according to the physics of photo- 
ionization presented in Section 15.11 and implemented it us- 
ing a single freq uency bin in th e parallel N-body-Tree-SPH 
code GADGET- 2 (|Springelll2005l '). The description of this im- 
plementation is the subject of this section. 



5.2.1 Transport of ionizing photons 

The transport of ionizing photons is performed in finite ra- 
diative transfer time steps of size At r , employing the absorp- 
tion coefficient k„ given by Eqs. 1141 and 1151 Photon packets 
are transported in cones according to the description given 
in Section 3] At the end of each radiative transfer time step, 
i.e. at time t r + At r , where t r is the current simulation time, 
we know the number of ionizing photons AMin,% impinging 
on and the number of ionizing photons AA^ts,, absorbed by 
particle i over the time interval At r (cp. Section [4.3|> . The 
photo-ionization rate Ti is obtained directly, that is without 
explicitly referring to the mean intensity J„, using 



m H 



*—YY At r = AN in ,i [1 - exp(-Tf-) , 



where Xi is the hydrogen mass fraction and 

ANabsA 



= - In 1 - 



AAfi„A 



(21) 



(22) 



is the a posteriori optical depth and we use superscripts to 
indicate the time at which quantities are evaluated. In the 
next section we describe how the photo-ionization rate is 
used to update the neutral fraction of particle i. 



5.2.2 Solving the rate equation 

In the simulation the neutral fraction associated with any 
gas particle i is assumed to be constant over the time step 
At r (see Section [4. 4p . The correspondingly discretized rate 
equation would read (cp. Eq. I17[l . 



t r + At r 



n, 



~At r 



(23) 



According to Eq. 1191 however, the neutral fraction evolves 
on the time-scale T eq during the photo-ionization process. 
In order to accurately follow the evolution of the neutral 
fraction using Eq. [23] one would have to choose time steps 
At r *C T eq limited by the ionization equilibrium time-scale, 
over which our assumption of a constant neutral fraction is 
a good approximation. 

In our implementation we choose the radiative trans- 
fer time step At r independently of the time-scale T eq , which 
can be prohibitively small, by employing the following sub- 
cycling strategy at the end of each radiative transfer time 
step. We only assume that the flux dMi n ,i/dt of ionizing 
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photons impinging on particle i is constant over the time 
step At r , i.e. cLN~% n ,i/dt = AJ\f in ,i/At r , in agreement with 
the temporal discretisation of the radiation transport in 
traphic. We then a posteriori follow the evolution of the 
neutral fraction in time ti £ (t r ,t r + At r ) on sub-cycles 
SU ^ At r , 



ti+Sti 



(24) 



Our assumption of a constant ionizing flux implies that the 
photo-ionization rate T oc (1— e~ T )rj~ 1 (see Eq.l21[). Hence, a 
change in the neutral fraction implies a change in the photo- 
ionization rate F 



r? = r 



1 - exp(-r ' 



exp( 



(25) 



where r* r and t\ t are the photo-ionization rate and the 
optical depth at the beginning of the sub-cycling, given by 
Eqs.[2l]and[2! and r* 1 = T^rfr /r)* r . 

The number of ionizations AjV 3U (,,i occurring over At r 
is then AN~ sub)i = {m t i r Xl r /m H ) J2 ^Vi^ti, where the 
sum is over all sub-steps Sti in (t r ,t r + At r ). We set 
8ti = min(/r** iy t r + At r — ti), where / < 1 is a dimen- 
sionless factor. It can be demonstrated that for this choice 
of Sti Eq. [24] is a sufficiently good approximation to Eq. [17] 
that the neutral fraction never violates the physical bound 
^ T)i ^ 1. Because the neutral fraction changes during the 
sub-cycling, the number of ionizations AN S ub,i can be less 
than the number of photons AN a bs,i that have been removed 
due to absorptions during the radiation transport over the 
time step At r based on the assumption of a constant neu- 
tral fraction. We then explicitly conserve photons by adding 
AJ\f a b s ,i — AjV s „5,, photons to the photon transport in the 
next radiative transfer step. 

When either the photo-ionization rate or recombination 
rate is high, r eq and hence St will be very small (dropping the 
particle index i for simplicity) . For St <C At r the sub-cycling 
becomes computationally very expensive. We find, however, 
that photo-ionization equilibrium is typically reached after 
only a few sub-cycles. Once photo-ionization equilibrium is 
reached, an explicit integration of the rate equation is no 
longer necessary. For a photon-conserving transport we still 
require knowledge about the number of photo-ionizations 
and recombinations occuring during the equilibrium phase. 
Both can, however, be obtained in a stroke, based on the 
number of photo-ionizations and recombinations occuring 
during the last sub-cycle step over which the rate equation 
was explicitly integrated. 

The importance of properly following the evolution of 
the photo-ionization rate in the presence of an evolving neu- 
tral fr action in the optically thiclf^l regime has been pointed 
out bv lMeilema et all (|2006T ) . There, a time-averaged photo- 
ionization rate obtained from an iterative procedure was em- 
ployed. The sub-cycling procedure (Eq. I25p presented here 
has the advantage that it is well-motivated also in the pres- 
ence of recombinations. 

We note that whenever one is only interested in obtain- 
ing the equilibrium neutral fraction, the detailed handling 
of the photo-ionization rate is rather unimportant. This is 

10 For r <g 1 Eq. 1251 implies that the photo-ionization rate is 
constant, T* s = r' r . 



because ionization equilibrium implies that the number of 
photo-ionizations cU\Ti„/dt(l — e~ T )At r over the time in- 
terval At r exactly cancels the number of recombinations 
(1 — rj) 2 NnriHOtAtr over that same time interval. This bal- 
ance, however, has a unique (and stable; see Eq. IT9|) solu- 
tion for the neutral fraction. The equilibrium neutral frac- 
tion then also implies the correct photo-ionization rate, via 
Eq. 1211 When one is interested in following the details of the 
evolution of the neutral fraction towards photo-ionization 
equilibrium, on the other hand, the dependence of the photo- 
ionization rate on the neutral fraction needs to be taken into 
account, as presented above. 



5.2.3 The time step At r 

Our main consideration when choosing the size of the radia- 
tive transfer time step for the simulations we are presenting 
in this work, is that we wish to accurately reproduce the 
analytical and numerical reference results we are comparing 
with. These results include the time-dependence of the size 
of ionized regions around ionizing sources. At early times, 
just after the sources start to emit ionizing photons, the 
ionized regions expand quickly into the neutral hydrogen 
field surrounding the sources. To accurately reproduce this 
early phase of rapid growth, we necessarily have to employ 
time steps At r that are relatively small. The phase of rapid 
growth is, however, only of relatively short duration. The 
subsequent evolutionary stages of modest resp. slow growth, 
which account for most of the (simulation) time, are often 
more interesting. We show in Section ["5.3.1l that whenever we 
are not interested in the very early phase of rapid growth 
we can, thanks to the photon-conserving nature of traphic, 
choose substantially larger time steps without affecting the 
outcome of our simulations. 

For all but one of the simulations we present in 
this work, we will be concerned with solving the time- 
independent radiative transfer equation. In these simula- 
tions, we choose to propagate photon packets downstream 
from their location of emission over only a single inter- 
particle distance per radiative transfer time step, unless 
stated otherwise. This approach is equivalent to solving the 
time- independent radiative transfer equation in the limit of 
small radiative transfer time steps, At r — » 0. We have ex- 
plicitly checked for all our simulations that the time step 
was chosen sufficiently small to be in agreement with this 
limit. 

Our treatment of the time-independent radiation trans- 
port reduces the computational effort for the simulation of 
problems for which the time-step has been fixed to a small 
value, e.g. by considerations like those presented in the be- 
ginning of this section. In the limit that radiation com- 
pletely fills the simulation box, the computational effort re- 
quired to solve the time-independent radiative transfer equa- 
tion over the time interval T by propagating photon pack- 
ets only over a single inter-particle distance per radiative 
transfer time step At r is proportional to (cp. Section |4. 2. 3|l 
Nsph x N ng b x N c x T/At r . This has to be compared to 
the computational effort required to follow all photon pack- 
ets over each time step until they leave the simulation box, 
which is proportional to Nsph x N ng b x 7V C x N^p H x T/At r . 
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This is larger by a factor of N S ' PH , which for typical simu- 
lations reaches values of the order of 100. 

In Sec. I5.3,3l we will present one simulation in which we 
solve the time-dependent radiative transfer equation. In this 
simulation we will employ the photon clock (Section 14, 4|l to 
accurately control the distance over which photon packets 
are propagated over each time step At r to match the light 
crossing distance cAt r . In this context it is interesting to 
note that in the case of small time steps, i.e. cAt r < Lt ox , it 
is less expensive to solve the time-dependent than the time- 
independent radiative transfer equation. This is because the 
computational effort to solve the time-dependent equation 
(again in the limit that radiation completely fills the sim- 
ulation box and assuming that its boundaries are photon- 
absorbing) scales as Ns phX N ngb xN c x cAt r /L box x Ngp H x 
T/At r , which is smaller than the computational effort for 
obtaining the time-independent solution (assuming that over 
each radiative transfer time step photon packets are trans- 
ported until they leave the box) by the factoJ"! cAt r /Lbox- 



5.2.4 Effective multi-frequency description 

In our current implementation we use only a single fre- 
quency. Thus, we either assume that the ionizing radi- 
ation is mono-chromatic, or we assume ionizing radia- 
tion with a frequency spectrum J v and provide an ef- 
fective multi-frequency description using only a single 
frequency bin (for a textbook introduction to the nu- 
merical treatment of multi-fre quency radiation, see e.g. 
iMihalas fc Weibel Mihalaslll984h . 

For the latter case, we define a frequency-independent 
(grey) photo-ionization cross-section a such that the total 
photo-ionization rate (Eq. I16[) is conserved, 



r = 



' 4jrJ„(i/) 
dv — - — -^(Jv = a 



dv 



a — dv — - — —t-Ou X 



dv 
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(26) 
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5.3 Tests 

In this section we report on the performance of our imple- 
mentation in simple and well-defined test problems. The 
tests comprise the evolution of the ionized region around 
a single star in a homogeneous medium (Section |5.3.ip . 
the casting of a shadow behind an opaque obstacle (Sec- 
tion l5.3.2|) and the propagation of ionization fronts driven by 
the ionizing radiation of multiple stars in an inhomogeneous 
density field obtained from a cosmological simulation (Sec- 
tion [5]3]3}. We compare the results obtained with traphic 
with analytic solutions, where available. For most physical 



11 Observe that for larger time steps, i.e. At r Lfj OX /c, and 
assuming photon-absorbing boundaries, the computational effort 
for solving the time-dependent radiative transfer equation equals 
the computational effort for solving the time-independent radia- 
tive transfer equation. It exceeds the computational effort for solv- 
ing the time-independent radiative transfer equation with photon 

packets propagating only a single intcr-particle distance per ra- 

1/3 

diative transfer time-step if cAt r > Lb ox /N gpjj- 



settings, however, no analytic solution to the radiative trans- 
fer problem is known, mainly due to the complex geome- 
tries involved. We have therefore designed the test problems 
in Sections 15.3.11 and 15.3.31 to closely follow the description 
given in the Cos mological Radiat ive Transfer Codes Com- 
parison Project (|lliev et al.ll2006l ). which provides a set of 
very useful numerical reference solutions to compare with. 

Throughout we will assume that the dens ity field is 
static . To facilitate a direct comparison with llliev et aU 
(2006), we present results after mapping physical quan- 
tities defined on the SPH particles to a regular grid, 
unless stated otherwise. This is done using a mass- 
co nserving SPH interpolation s i milar to the one described 
in lAlvarez. Bromm. fc Shapiro! (|2006I ). We opted for the 
SPH interpolation since it is most consistent with the 
SPH simulation we are employing. For comparison, we re- 
peate d our analysis using TSC m ass-conserving interpola- 
tion (jHocknev fc Eastwood! [19 88) but found no significant 
differences. 

For all tests reported in this section we em ploy the on- 
the-spot approximation (e.g. IOsterbrocklll989h. in order to 
allow a direct comparison with llliev et alj 1)20061) . In the 
on-the-spot approximation diffuse photons emitted during 
recombinations to the hydrogen ground energy level are as- 
sumed to be immediately re-absorbed by neutral hydrogen 
atoms close to the location of emission. The effect of dif- 
fuse recombination radiation can then be simply taken into 
account by setting the recombination coefficient a to the so- 
called case B value a_g, which for the temperature range of 
interest can be well approximated by 



qs(T) = 2.59 x 10" 
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(28) 



We will report on the study of diffuse radiation in which 
we will be explicitly following the scattering of diffuse pho- 
tons instead of employing the on-the-spot approximation in 

future work. 

In their simulations, llliev et alj (|2006l ) assumed that 
the speed of light is infinite, i.e. they solved the time- 
independent radiative transfer equation. For the comparison 
with these simulations we will therefore make the same as- 
sumption (recall Sec. l5.2.3l for the discussion of how we solve 
the time-independent radiative transfer equation). From 
now on, when referring to the radiative transfer equation, 
we therefore assume it to be of the time-independent form, 
unless stated otherwise. We will repeat one of the simula- 
tions presented in Sec. 15.3.31 to solve the time-dependent 
radiative transfer equation by employing the photon packet 
clocks as described in Sec. 14.41 

5.3.1 Test 1: Spherically-symmetric HII region expansion 

In this section we consider the problem of a steady ion- 
izing source emitting My mono-chromatic photons of fre- 
quency h p v — 13.6 eV per unit time, which is turned on 
in a static, initially neutral, uniform medium of constant 
hydrogen number density tih- 

The analytical solution to this p roblem is well-known 
(for a textbook discussion see e.g. iDopita fe Sutherland! 
2003). In equilibrium, the number of photons emitted by 
the source has to compensate the number of recombina- 
tions within the surrounding, stationary, ionized region, the 
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Figure 6. Test 1. Slice through the simulation box at z = Li, ox /2 showing the neutral fraction at the end of the simulation (t r = 500 Myr), 
i.e. in photo-ionization equilibrium, for simulations with (bottom row) and without (top row) resampling of the density field. The angular 
resolution increases from N c = 8 in the left-most to N c = 128 in the right-most column, as indicated in the panel titles. The spatial 
resolution is fixed to Nsph = 64 3 , N ng b = 32 and is indicated by the cross of length 2(h) in the upper left corner of each panel. Black 
contours show neutral fractions of r\ = 0.9, 0.5, logr; = —1, —1.5, —2, —2.5, —3, —3.5, —4, going from the outside in. The white dot-dashed 
circle indicates the Stromgren sphere, which should be, and is, just inside to the r) = 0.5 contour. The colour scale is logarithmic and has 
a lower cut-off log ?; = —5. Note that the resampling strongly suppresses the particle noise seen in the top-row panels. 




Figure 7. Test 1. Spherically averaged neutral and ionized fractions within the Stromgren sphere at t r = 500 Myr for different angular 
resolutions, as indicated in the legend. The profiles in the left-hand (right-hand) panel correspond to simulations without (with) resampling 
the density field. The spatial resolution is fixed (Nsph = 64 3 , N ng b = 32). The black solid curves corresponds to the exact pr ofiles, given 
by the solution of Eq. 1331 The grey bands show the range of neutral and ionized fractions found by other codes as reported in llliev et al.l 
(2006). The vertical dotted line marks the radius r = 2(h), corresponding to the spatial resolution employed. In the right-hand panel, 
the additional (right-most) vertical dotted line indicates the radius corresponding to the spatial resolution 2(h) of the SPH simulations, 
which is the scale on which particle positions are randomly displaced during the resampling. 
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Figure 5. Test 1. The evolution of the ionization front for the 
angular resolutions N c = 8, 16, 32, 64 and 128, as indicated in the 
legend. The spatial resolution is fixed (Ngpn = 64 3 , N ng b = 32). 
The top panel shows the position of the ionization front r/, nIlm 
normalized by the Stromgren radius r a as a function of time. The 
thick black solid curve shows a numerical reference solution ob- 
tained with a one-dimensional, grid-based radiative transfer code 
(see text). The grey curve shows the analytic reference solution, 
Eq. 1321 which has been obtained by assuming \ = 1 throughout 
the ionized region. The results from the numerical simulations em- 
ploying TRAPHIC closely match the numerical reference solution. 
The bottom panel shows the position of the ionization fronts of 
the top panel divided by the analytic reference solution. Note that 
the analytic reference solution slightly differs from the numerical 
reference solution, due to the simplifying assumptions inherent to 
the analytic approach (see also the discussion of Eq. I33( . 



so-called Stromgren sphere. Assuming that the Stromgren 
sphere is fully ionized, i.e. \ = lj its radius r s is therefore 
given by the balance equation 



■My = ^r z a oi B {T)n 2 H . 



(29) 



The evolution of the spherical, fully ionized region to- 
wards the equilibrium Stromgren sphere is described by the 
following equation for its radius rj, the ionization front, 



a 2 dri \r 4 3 fr\ 2 

AnrjUH—jj- = My — -nr I aB{-i )n H . 



(30) 



Introducing the new variables £ = ri/r s and r = t/r s , where 
the Stromgren time-scale t s = l/^aBnu) is the recombina- 
tion time in a fully ionized gas, we arrive at the differential 
equation 

(31) 



dr 



3£ 2 



the solution of which reads 



n(t) = r 3 (l 



t/T 3 S.l/3 



(32) 



Hence, the ionization front reaches the Stromgren sphere 
after a few Stromgren times r s and stays static thereafter. 

In reality the neutral fraction inside the ionized region 
is, however, not zero, but varies smoothly with the distance 
to the ionizing source. We therefore invoke the commonly 
employed definition of the ionization front as the radius at 
which the neutral fraction equals 50 per cent, i.e. rj — 0.5. 
The equilibrium neutral fraction rj eq (r) = lim^oo v(r) can 
be obtained by solving the equation (e.g. IOsterbrocklll989l ) 



Tjeq{r)nH 

4irr 2 



du A/" 7 (^)e T "o y = Xeq( r ) n H a B, 



(33) 



(34) 



where the optical depth T„{r) is given by 

r„(r) = n H (j v I dr 1 ?7 e «j(r'). 
Jo 

We refer to the neutral fraction rj eq (r) obtained by direct 
numerical integration of Eq. [33] as the exact neutral frac- 
tion profile and to Xeg{r) = 1 — ileq{ r ) as tne exact ion- 
ized fraction profile. These profiles are shown as black solid 
curves in Fig. [7] Note that for our choice of parameters, 
r 'i,eq = 1-05 r s . The ionization front obtained from the solu- 
tion to Eq.[53]is thus at a slightly larger radius than the equi- 
librium ionization front obtained assuming the Stromgren 
sphere is fully ionized (Eq. I32|l . 

The last observation indicates that the analytic solu- 
tion given by Eq. [32] generally fails to provide an accurate 
reference solution that can be used to judge the validity of 
the numerical results obtained with a new radiative trans- 
fer scheme like TRAPHIC, due to its simplification of the 
problem. We therefore employ a one-dimensional, explic- 
itly photon-conserving, grid-based radiative transfer scheme 
that we have implemented ourselves to obtain an accurate 
numerical reference solution. We will make use of this nu- 
merical reference solution in our comparisons below. 

In the following we present a suite of radiative trans- 
fer simulations demonstrating that TRAPHIC is able to accu- 
rately follow the evolution of the ionization front around a 
single ionizing point source. For the setup of the numer- 
ical simulations we closely follow the description of Test 
1 presented in llliev et ah (|2006h . the only difference be- 
ing the spatial resolution we employ. This allows us to di- 
rectly compare our results to the results presented in the 
code comparison project. In particular, we choose a num- 
ber density «h = 10~ 3 cm" 3 and an ionizing luminosity 
of A/" 7 = 5 x 10 48 photons s _1 . The gas is assumed to 
have a constant temperature T = 10 4 K. With these val- 
ues, r s — 5.4 kpc and r s — 122.4 Myr. 

We set up the initial conditions in a simulation box of 
length Lbox = 13.2 kpc containing Nsph = 64 3 SPH parti- 
clea 12 !. The ionizing source is located in the centre. The box 
boundary is photon-transmissive. We assign each SPH par- 
ticle a mass m = numu L Box /Nsph ■ The positions of the 
SPH particles are chosen to be glass-likf 3 ! Glass-like initial 
conditions yield a more regular distribution of the particles 
within the box as compared to Monte Carlo sampling of 
the density field. The SPH smoothing kernel is computed 
and the SPH densities are found using the SPH formalism 
implemented in GADGET-2, with N ng b = 48. 

We perform 5 simulations, increasing the angular res- 
olution in factors of two from iV c = 8 to N c = 128. The 
number of neighbours employed for the transport of radi- 
ation is fixed to N ng t = 32. Hence all 5 simulations em- 
ploy the same spatial resolution. The time step is set to 

12 We note that llliev et al.l ll2006h employed N cM = 128 3 cells, 
with the ionizing source located in one of the corners of the box. 

13 The glass-like distribution of particles is achieved by first plac- 
ing them randomly in the simulation box and thereafter letting 
them freely evolve under the influence of a reversed-sign (i.e. re- 
pelling) gravitational f orce unt i l they settle down into an equilib- 
rium configuration, see 
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At r = 10 4 yr. In Fig. [5] we show the evolution of the ion- 
ization front radius, which we determined by taking the av- 
erage over the positions of all particles that have a neutral 
fraction 0.4 < rj < 0.6. For all 5 simulations, the position 
of the ionization front never deviates more than 5 per cent 
from the analytic solution, Eq. 1321 comparable to what has 
been found with other codes as reported in the Cosmologi- 
cal R adiative Transfer Code Comparison Project (jlliev et al.l 
2006). Note that the deviations from the analytic solution 
can mainly be attributed to the fact that the analytic ap- 
proach provides only an approximate expression for the ra- 
dius of the ionization front, because it erroneously assumes 
X = 1. In fact, all simulations (except for the N c = N ng t run, 
which shows a small deviation of less than two percent, the 
reason for which will become clear in our discussion below) , 
nearly perfectly follow the numerical reference solution and 
approach the proper asymptotic limit rr,eg = 1-05 r s . 

The top row of Fig. [6] shows the neutral fraction in 
a slice through the centre of the simulation box at t T = 
500 Myr, which marks the end of the simulation, for each 
of the 5 simulations. As we already noted, the ionization 
front is at the expected position. As is true for all other 
surfaces of constant neutral fraction shown, the ionization 
front clearly exhibits the expected spherically symmetric 
shape, although it is noisy in some of the simulations. The 
amount of noise depends on the ratio of the angular and 
spatial resolution employed. For N c = 8 (left-most panel in 
the top row), the average number of neighbours per emis- 
sion or transmission cone is high, N ng b/N c = 4 and, as a 
result, numerical noise arising from the representation of 
the continuous density field with discrete SPH particles are 
suppressed. With increasing angular resolution the average 
number of neighbours per cone decreases, and the contours 
become more noisy. The noise level reaches a maximum at 
JV C = N ng b (middle panel in the top row). For higher an- 
gular resolutions, the probability of finding no neighbours 
inside an emission or transmission cone becomes high and 
a large number of ViPs are created. The ratio of the num- 
ber of ViPs to the number of SPH particles enclosed by the 
ionization front for the simulation with angular resolution 
N c = 8,16,32,64 and 128 is » 0,0.003,0.06,0.5 and 0.9, 
resp. The ViPs placed in empty cones regularize the neutral 
fraction of the ionized density field by distributing the pho- 
tons they absorb amongst their N ng b SPH neighbours using 
(photon-conserving) SPH interpolation. 

In the left-hand panel of Fig. [7] we plot the neutral and 
ionized fraction averaged in spherical shells as a function of 
distance to the star, for all 5 simulations. The grey bands 
indicate the neutral and ionized fraction profiles obtained 
with other codes as reported in t he cosmological radiative 
transfer code comparison project (|lliev et al.ll2006l ). Except 
for the N c = 32 run, for which the neutral contours were 
most noisy (see Fig. [6} , all simulations agree very well with 
the results published in the comparison project. The devi- 
ations from the exact equilibrium neutral fraction profile, 
i.e. the result of the numerical integration of Eq. 1331 can be 
explained by the particle noise seen in Fig. [S] Due to the 
fuzzy structure exhibited by the neutral fraction contours, 
a range of neutral fractions can simultaneously be found 
within each spherical shell. The profiles obtained from the 
numerical simulation with traphic should therefore not be 
directly compared to the exact profile, i.e. the solution of 



Eq. 1331 but to the profile that results after locally averaging 
the exact profile along the radial direction. 

To investigate the effect of particle noise on the neu- 
tral fraction profile we apply the resampling technique in- 
troduced in Section ^. 51 We perform a series of 5 simulations 
that are identical to the simulations presented above, except 
that every 10th radiative transfer time step the particle po- 
sitions are randomly perturbed within their SPH spheres of 
influence. The numerical implementation of the resampling 
is described in Appendix IB3I The densities are not recal- 
culated after the positions have been changed due to the 
resampling, because this would generate fluctuations in the 
neutral hydrogen density which would increase the recombi- 
nation rate and lead to a smaller Stromgren sphere. The re- 
sulting neutral fraction profiles are shown in the right-hand 
panel of Fig. [7] All profiles are now in close agreement with 
each other and with the exact result. The panels in the bot- 
tom row of Fig.|S]show the neutral fraction in a slice through 
the centre of the simulation box from the simulations with 
resampling. Clearly, resampling strongly suppresses the par- 
ticle noise visible in the panels in the top row of Fig. [5] 
yielding nearly spherical neutral fraction contours. 

We now investigate the dependence of the equilibrium 
neutral and ionized fraction profile on the spatial resolu- 
tion by varying Nsph, the number of particles employed 
in the simulation. Because traphic is explicitly photon- 
conserving, we expect that the radiative transfer in a ho- 
mogeneous medium is essentially inde pendent of the spatia l 
resolution (see e.g. the discussion in iMellema et all 120061 ) . 
except for a trivial averaging. For each of the simulations 
with angular resolution N c = 8, 32 and 128 and N ng b = 32 
presented above, we performed three additional simulations, 
at lower (N S ph = 16 3 ,32 3 ) and higher (N S ph = 128 3 ) 
spatial resolutions, but otherwise identical to the fiducial 
(Nsph = 64 3 ) case. We will focus on the N c = 8 runs, but 
note that the N c = 32 and N c = 128 series shows the same 
behaviour. 

The equilibrium neutral and ionized fraction profiles are 
shown in the left-hand panel of Fig. [5] For all spatial reso- 
lutions the ionization front is at nearly the same radius. It 
can furthermore be seen that when the spatial resolution is 
increased, the equilibrium neutral fraction follows the exact 
result more closely. The simulation employing the fiducial 
spatial resolution (Nsph = 64 3 ) is almost converged. The 
differences in the neutral fraction profiles between the simu- 
lations using different numbers of particles are fully consis- 
tent with the corresponding spatial resolutions, as is demon- 
strated in the right-hand panel Fig. [5] There, we average 
the neutral fraction profiles obtained in the simulations em- 
ploying Nsph = 32 3 ,64 3 and 128 3 particles over the spa- 
tial resolution element of the lowest resolution simulation 
(Nsph = 16 3 ), the size of which is indicated by the vertical 
line. The averaged profiles match the neutral fraction pro- 
file obtained in the low spatial resolution simulation almost 
exactly. This shows that decreasing the spatial resolution 
does not introduce any artefacts. The solution obtained by 
traphic is the converged solution averaged over the adopted 
spatial resolution. 

Finally, we show how the size of the time step At r 
affects the outcome of our simulations. We again concen- 
trate on the simulation with angular resolution 7V C = 8 
(and Nsph = 64 3 , N ng b = 32), noting that the simulations 
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Figure 8. Test 1. Spherically averaged neutral and ionized fractions within the Stromgren sphere at t r = 500 Myr. The simulations 
all have the same angular resolution (N c = 8) and employ the same number of neighbours (Nngb = 32), but use a different number of 
SPH particles, increasing from N SPH = 16 3 to N SPH = 128 3 in factors of 2 3 . The black solid curves corresponds to the ex act profiles, 
obtain ed from Eq. 1331 The grey bands show the range of neutral and ionized fractions found by other codes as reported in llliev et al.1 
(2006). Left-hand panel: For each simulation the spatial resolution is marked by a vertical line (in colour and line-style identical to the 
corresponding profile) at radius 2(h). The higher the spatial resolution, the more closely the exact profile is approached. Right-hand 
panel: The profile corresponding to the lowest spatial resolution simulation (Ngpu = 16 3 , blue dot-dashed line) is repeated from the 
left-hand panel. The profiles of all other simulations have been averaged over the spatial resolution element 2(h) of the lowest spatial 
resolution simulation, corresponding to the radius marked by the vertical line. The profiles become nearly identical after smoothing them 
to the same resolution. 
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Figure 9. Test 1. Effect of the size of the time step. We show the 
evolution of the ionization fronts for four simulations with N c = 8, 
N SPH = 64 3 , N ngb = 32 and time steps At r = 0.01,0.1, 1 and 
10 Myr, resp, as indicated in the legend. After an initial phase, 
the evolution of the ionization fronts becomes independent of the 
size of the radiative transfer time step. 

of higher angular resolution exhibit a similar behaviour. In 
Fig. [9] we show the evolution of the ionization front for four 
different choices for the size of the radiative transfer time 
step, At r = 0.01,0.1,1 and 10 Myr. Note that the simu- 
lation with At r = 0.01 is identical to the simulation dis- 



cussed in Fig. [5] In order to keep the angular sampling the 
same, at each radiative transfer step we split the emission 
of photons over 10, 100 and 1000 random orientations of the 
emission cone tessellation of the ionizing source for the sim- 
ulations employing At r = 0.1, 1 and 10 Myr, resp. Photon 
packets that are emitted by the source in a certain orienta- 
tion are transmitted further downstream and can propagate 
over multiple inter-particle distances. We follow each photon 
packet until it has either been absorbed or left the simula- 
tion box, to properly solve the time-independent radiative 
transfer equation for the large time steps under considera- 
tion. 

From Fig. [5] we see that the evolution of the ionization 
front for the simulations with time step At r = 0.1, 1 and 
10 Myr is delayed with respect to the evolution of the ioniza- 
tion front for the simulation with time step At r = 0.01 Myr. 
This delay increases with the size of the time step, being 
barely visible for the simulation using At r = 0.1 Myr. The 
delay arises because the neutral fraction is only updated at 
the end of each radiative transfer time step. Photons that 
have been absorbed during the transport over a single time 
step but that have not been used to advance the neutral 
fraction during the subsequent sub-cycling of the rate equa- 
tion are therefore only re-inserted in the transport process 
at the beginning of the next time stetPL and are thus de- 



We note that we have also tried to re-insert these photons 
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layed. From Fig. [9] it can, however, be seen that after a few 
time steps (here w 15), the ionization front catches up to 
agree with the ionization front obtained in the simulation 
using At r = 0.01 Myr. We have convinced ourselves that 
from then on the neutral fraction profile around the ioniz- 
ing source is also nearly identical to the profile obtained in 
the simulation with At r = 0.01 Myr (Fig. [7j). 

In summary, in this section we showed that traphic 
is able to reproduce the expected equilibrium neutral frac- 
tion around an ionizing source embedded in a homogeneous 
medium, as well as the dynamical evolution towards it. Be- 
cause the radiative transfer is explicitly photon-conserving, 
the spatial resolution only determines the scale over which 
the converged solution is averaged. We have furthermore 
seen that the performance of traphic is stable with respect 
to variations in the size of the time step. Particle noise due 
to the discrete nature of SPH simulations is small except for 
the choice of parameters N c = N ng b- The noise can be suc- 
cessfully suppressed by applying the resampling technique 
of Section [4.51 i.e. by periodically perturbing the positions 
of the SPH particles within their spheres of influence. 

In the next section we employ a test with broken sym- 
metry to demonstrate the ability of traphic to correctly 
produce shadows behind opaque obstacles and we study fur- 
ther how the choice of the parameters N c and N ng b affects 
the outcome of the numerical simulations. 

5.3.2 Test 2: Breaking the spherical symmetry 

In the absence of scattering interactions, photons propagate 
along straight lines into the direction set at the time of their 
emission. Consequently, opaque obstacles throw sharply de- 
fined shadows. In this section we are mainly interested in 
studying the properties of the shadow thrown by an opaque 
obstacle exposed to ionizing radiation from a single point 
source, as obtained with traphic. At the same time, we will 
extend the study of particle noise presented in Section r5.3.1l 
to include other choices for the parameter N ng b. 

The geometry of our test p roblem closely follow s the 
description of the shadow test in iMellema et all (|2006T ) . We 
consider a source emitting 7v 7 = 10 54 photons s" 1 , each of 
hydrogen-ionizing energy h v v = 13.6 eV, residing in an ini- 
tially neutral, static hydrogen-only field of constant number 
density riu = 1.87 x 10 -4 cm -3 and temperature T = 10 4 K. 
The Stromgren radius (Eq. I29J1 corresponding to this set of 
parameters is r s = 0.965 Mpc and the Stromgren time is 
T a = 654.3 Myr. For the numerical simulation a star particle 
is placed in the centre of a box of size Lbox = 1 Mpc. The 
boundaries of the box are photon-transmissive. The density 
field is sampled using Nsph = 64 3 gas particles with mass 
m = nHrriHL% ox /NspH at glass- like positions. The parti- 
cle densities are found using the SPH interpolation imple- 
mented in GADGET-2, with N ngb = 48. We furthermore place 

immediately after they have been absorbed, by integrating the 
rate equation already at the end of each transport cycle (without 
updating the neutral fraction) to obtain the number of photons 
that have been erroneously counted for being absorbed (because 
of the assumption of a constant neutral fraction). The re-insertion 
of these photons in the transport process within the same radia- 
tive transfer time step over which they have been absorbed indeed 
reduced the delay of the ionization front observed in Fig. [9] 



an infinitely thin opaque disc perpendicular to the rc-axis at 
a distance of 0.08 Mpc along the a;-axis from the star (thick 
blue vertical lines in Figs. 1 1UI 14[) . The y and z coordinates 
of the disc centre are identical to those coordinates of the 
star. Photons that cross the disc are removed. 

We performed a series of radiative transfer simulations 
(with time step At r = 10 4 yr), using different choices for 
the parameters N ng b, which sets the spatial resolution if the 
total number of SPH particles is fixed, and N c , which sets 
the angular resolution. The results are shown in Fig. 1101 - 1141 
displaying a slice through the simulation box at z = Lbox/2- 
In each panel, the black dash-dotted line shows the expected 
position of the ionization front (Eq. l32p at time t r — 80 Myr, 
which marks the end of the simulation. The black solid lines 
emerging from the star at the centre of the slice show the 
boundaries of the shadow expected to be thrown by the 
opaque disc (thick blue vertical line). In the top-left corner 
of each panel we indicate the spatial resolution by a cross 
of length 2{h) corresponding to the average diameter of the 
sphere containing N ng t neighbours. The angular resolution 
is indicated by the angle u) enclosed by the black dashed line 
and the upper shadow boundary, where uj is determined us- 
ing Eq. [9] It indicates the maximum angle photons can theo- 
retically diverge from the shadow boundary into the shadow 
region, given the chosen angular resolution N c . 

The position of the ionization front (the iso-surface for 
which the neutral fraction r\ = 0.5) at time t r = 80 Myr 
as obtained with traphic is shown by the red solid line. In 
all panels, that is for all spatial and angular resolutions, the 
ionization front is found at the proper position, in agreement 
with our findings of Section r5.3.1l The shadow thrown by the 
opaque disc is always sharp. We now discuss the dependence 
of the results on the chosen spatial and angular resolutions. 

In Fig.QJj]we show the ionization front obtained in sim- 
ulations employing a fixed spatial resolution, N ng b = 32, but 
varying angular resolution, ranging from N c = 8 in the left- 
most to N c = 128 in the right-most panel. The most promi- 
nent difference between the results of the different simula- 
tions is the noisiness of the contour tracing out the ionization 
front. The angular resolution study presented here is very 
similar to the one in the last section. For the lowest angular 
resolution, N c = 8, the ionization front is very smooth due 
to the large number of neighbours within each emission and 
transmission cone. The noise increases with the angular reso- 
lution until it reaches a maximum for N c — N ng b. For higher 
angular resolutions, particle noise is efficiently suppressed 
due to the large number of ViPs that are placed in empty 
cones and that distribute the photons they absorb amongst 
their N ng b SPH neighbours using (photon-conserving) SPH 
interpolation. 

From Fig. [10] it can furthermore be seen how the sharp- 
ness of the shadow thrown by the opaque disc depends on 
the angular resolution. For the lowest angular resolution, the 
shadow is somewhat blurred, though not nearly as much as 
the angular resolution would imply. Increasing the angular 
resolution yields slightly sharper shadows. However, if the 
angular resolution is increased beyond N c — N ng b, the shad- 
ows become slightly less sharp. This is because the photons 
absorbed by ViPs are distributed amongst the neighbour- 
ing gas particles using SPH interpolation and the interpola- 
tion procedure does not know about the shadow region. The 
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Figure 10. Test 2. Slice through the simulation box at z = L box /2 showing the ionization front (red solid line) at time t r = 80 Myr 
around an ionizing source sitting in the centre of the simulation box. The black dot-dashed line shows the expected ionization front 
position. The thick blue vertical line indicates an obstacle opaque to ionizing photons and the black solid lines trace out the boundaries 
of the shadow this obstacle is expected to throw. The cross and the black dashed line indicate the spatial and angular resolution, 
respectively, as described in the text. The spatial resolution is fixed to Nsph = 64 3 , N ng b = 32. The angular resolution increases from 
N c = 8 in the left-most panel to N c = 128 in the right-most panel, in factors of 2. 
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Figure 11. Test 2. Same as Fig. 1101 but with the angular resolution fixed to N c = 32 and the spatial resolution decreasing, from 
Nsph = 64 3 , N ng b = 8 in the left- most panel to Ngpfj = 64 3 , N ng j, = 128 in the right-most panel, increasing the number of neighbours 
N ng b in factors of 2. 



slight diffusion of photons across the shadow boundary is in 
this case consistent with the spatial resolution. 

In Fig. [11] we show the results of the simulations where 
we fixed the angular resolution to N c =32, but varied the 
spatial resolution by changing N ng b- The trends visible in 
Fig. [10] can also be observed here. The ionization front is 
most noisy and the shadow is sharpest for N c = N ng t- For 
N ng b > N c , noise due to the discreteness of the particles 
employed for the transport of photons is suppressed by the 
large number of neighbours per cone, but the shadow is 
slightly blurred. The shadow becomes sharper for a smaller 
number of neighbours, since generally not all of the solid 
angle will be covered by the neighbours, an effect that be- 
comes more important for smaller numbers of neighbours. 
For N ng t < N c the ionization front becomes smoother due 
to the regularizing effect of the SPH interpolation from ViPs, 
which also leads to a small diffusion of photons across the 
shadow boundary, consistent with the spatial resolution. 

In Figs. [12] and [13] we keep the ratio N ng f, /N c fixed at 
N ng b/N c — 2 and N ng b/N c = 1/2, resp. In the first case 
there are on average 2 neighbours per cone, whereas in the 
second case there is on average one neighbour in every sec- 
ond cone. From Fig. [12] it is clear that the shadow does get 
sharper when the angular resolution is increased, although 
the effect is small, since the shadow is always very sharp. 
Because we keep the ratio N ngb /N c fixed at 2, the number 



of ViPs employed in the simulation stays low for all angu- 
lar resolutions and the shadows are not visibly diffused by 
the SPH interpolation of absorbed photons from the ViPs. 
Furthermore, the noisiness of the ionization front remains 
constant throughout the parameter range. This is because 
the noise is primarily set by the ratio N ng b/N c if N ng i, > N c . 
In Fig. 1131 on the other hand, there is a substantial proba- 
bility for creating a ViP per cone. Since the absolute number 
of ViPs present in the simulation increases with increasing 
angular resolution N c , the noise decreases with N c . 

In the last section (see bottom row of Fig. [(J} we em- 
ployed a resampling technique to decrease the noise exhib- 
ited by the neutral fraction contours. Recall from Section 
I4.2.2l that the apexes of the transmission cones are attached 
to the positions of the SPH particles. Hence, resampling 
results in slight shifts in the position of each transmission 
cone apex, on the scale of the spatial resolution employed 
in the SPH simulation. This shifting is expected to lead 
to a small diffusion of photons across the expected shadow 
boundary. As noted earlier in Section 14.51 such a diffusion 
due to particle resp. apex motion will also occur in radiation- 
hydrodynamical simulations because of the movement of the 
SPH particles. It is therefore interesting to study the prop- 
erties of the shadow thrown by an opaque obstacle in the 
case of resampling. 

In Fig. [I4]we show the results of a simulation which em- 
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Figure 12. Test 2. Same as Fig. 1101 but fixing the ratio between spatial and angular resolution to N ngb /N c = 2 (TVg 



■■ 64 3 ). 
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Figure 13. Test 2. Same as Fig. 1101 but fixing the ratio between spatial and angular resolution to N ng b/N c = 1/2 (Nsph 
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Figure 14. Test 2. Left-hand panel: Identical to the middle panel 
of Fig. HOI and llll Right-hand panel: Same as left-hand panel, ex- 
cept for the fact that in the simulation we periodically resampled 
the density field, resulting in a strong suppression of the particle 
noise seen in the left-hand panel. Note the (small amount of) dif- 
fusion of photons across the shadow boundary due to the motion 
of the transmission cone apexes. 



ploys the same parameters as used for the simulation pre- 
sented in the middle panels of Fig. \W\ and 1111 but with 
an additional resampling of the hydrogen density field ev- 
ery 10th radiative transfer time step. As in the last section, 



the resampling is performed without changing the hydrogen 
density, since this would lead to an enhanced recombination 
rate. Numerical noise is successfully suppressed by the ran- 
dom perturbations given to the positions of the SPH parti- 
cles. Consequently, the ionization front appears significantly 
smoother. The degree of photon diffusion into the shadow 
region is small and does not significantly degrade the angu- 
lar resolution of the radiative transfer. This is because the 
diffusion scale is set by the spatial resolution employed in 
the SPH simulation. Therefore, well-defined shadows will be 
thrown as long as the obstacle is spatially resolved. 

In summary, in this section we showed that traphic is 
able to produce a well-defined shadow behind an opaque ob- 
stacle, with the shadow sharpness in full agreement with the 
chosen spatial and angular resolutions. In fact, the shadows 
are much sharper than implied by the formal angular resolu- 
tion, thanks to the angular adaptivity inherent to traphic. 
For a fixed angular resolution, the shadows are sharpest for 
N c = N n gb- They are slightly broadened by photon diffu- 
sion for both N c < N ng b and N c > N ng b, due to the in- 
creased coverage of the solid angles traced out by the trans- 
mission cones with SPH particles for an increasing number 
of neighbours N ng b and the SPH interpolation of the pho- 
tons absorbed by ViPs, resp. We confirmed our finding of 
the last section that unless N ng b = N c , noise due to the dis- 
creteness of the particles on which the transport of photons 
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takes place is small, since it is suppressed by either the large 
number of neighbours per cone (if N c < N ng b) or the large 
number of ViPs employed (if N c > N ng b). The resampling 
technique presented in Section 14.51 is very effective at sup- 
pressing particle noise. We have seen that resampling the 
density field does not severely degrade the angular resolu- 
tion, even though it leads to a small shift of the cone apexes. 
As long as the opaque obstacle is spatially resolved by the 
SPH simulation, a well-defined shadow will still be thrown. 

The ability to produce sharp shadows is one of the main 
requirements a radiative transfer code has to pass. The re- 
sults of this section, together with the results of Test 1, 
which showed that traphic is able to reproduce the ex- 
pected neutral fraction within a spherically symmetric HII 
region, indicate that traphic can be used to perform the 
transport of ionizing photons in arbitrary complex geome- 
tries. This is the subject of Test 3, which we will describe 
next. 



5.3.3 Test 3: Cosmological density field 

In this test we study the propagation of ionization fronts 
around multiple sources in a static cosmological density 
field. The parameters of this test are taken from Test 4 
of the C osmological Rad iative Transfer Code Comparison 
Project (|lliev et al.| [2006'l. For reference we repeat the test 
description here. 

The initial conditions are provided by a snapshot (at 
redshift z « 8.85) from a cosmological N-body and gas- 
dynamical simulation performed using the c osmo logical 
(uniform-mesh) PM+TVD code of iRvu et all ljl993ft . The 
simulation box is L\, ox = 0.5 h~ x comoving Mpc on a side, 
uniformly divided into N ce ll = 128 3 cells. The initial tem- 
perature is fixed at T — 100 K everywhere. The halos in the 
simulation box were found using a friends-of-friends halo 
finder with a linking length of 0.25. The ionizing sources 
are chosen to correspond to the 16 most massive halos in 
the box. We assume that these have a black-body spectrum 
B v (u,T) with temperature T = 10 5 K. The ionizing pho- 
ton production rate is assumed to be constant and assigned 
assuming that each source lives for t a = 3 Myr and emits 
/ 7 = 250 photons per atom during its lifetime. Hence, the 
number of ionizing photons emitted per unit time is 

iloTTlBts 

where M is the total halo mass, Slo = 0.27, Qb = 0.043 and 
h — 0.7. For simplicity, all sources are assumed to switch 
on at the same time. The boundary conditions are photon- 
transmissive. Outputs are produced at t = 0.05,0.1,0.2,0.3 
and 0.4 Myr. 

With respect to the original test setup described above, 
we require three changes. First, since our code does not yet 
solve for the temperature of the gas, we assume a constant 
temperature of T = 10 4 K for the ionized gas. Second, since 
our code currently treats only a single frequency (bin), we 
assume the grey photo-ionization cross-section Eq. 1271 for 
which we find o — 1.49 x 10 -18 cm~ 2 . The third change con- 
cerns the input density field. Since our code works directly 
on the set of particles used in SPH simulations, we have to 
Monte Carlo sample the original input density field in order 
to place particles in the box. We replace every grid cell i by 




Figure 15. Test 3: Slice through the density field at z = Lf, ox /2. 
Contours show a neutral fraction of n = 0.5 (left-hand column) 
and r) = 0.01 (right-hand column) at times t = 0.05 Myr (top row) 
and t = 0.2 Myr (bottom row). Red contours show the results of 
our fiducial (N c = 32, N ngb = 32) simulation. For comparison, we 
sh ow the results of C 2 -RAY (green) and CRASH (blue), as reported 
in llliev et al. The agreement is excellent. 



Nsph — Mi/m SPH particles (randomly distributed within 
the volume of the grid cell), where Mi = piLl ox /N ce u is the 
mass of the cell and m is the mass of an SPH particle. If 
Ng PH is not an integer, we draw a random number from 
a uniform distribution on the interval (0,1) and place an 
additional particle if this number is smaller than the differ- 
ence between Ng PH and the nearest lower integer. We use 
Nsph = N ce u = 128 3 . Since the Monte Carlo sampling only 
results in the approximate equality N % SPH w Nsph, wc 
adjust the particle masses a posteriori to conserve mass, i.e. 
TTi — » 171 x Nsph / N SPH . After the particles have been 
placed, we calculate their densities using the SPH formalism 
of GADGET-2, with N ngb = 48. 

Note that Monte Carlo sampling the density field with 
Nsph — N ce u particles yields a smaller effective resolution 
than that of the grid input field in low density regions (many 
grid cells will be left empty of particles), and to a spurious 
higher resolution in high density regions (cells are sampled 
with many particles, even though there is no substructure on 
the scale of a single cell in the input field). Note also that 
because the initial conditions were specified on a uniform 
grid, we do not benefit from the intrinsic spatial adaptivity 
of TRAPHIC, effectively wasting computational resources. 

We performed a set of three radiative transfer simula- 
tions with angular resolutions N c = 8, 32 and 128, which 
we refer to as the low angular resolution, fiducial and high 
angular resolution simulations, resp. Every simulation used 
N ngb = 32 neighbours. The time step was set to At r — 
100 yr. 

In Fig. [K]we show in red the neutral fraction contours 
rj = 0.5 (left-hand column) and rj = 0.01 (right-hand col- 
umn) at times t = 0.05 Myr (top row) and t — 0.2 Myr (bot- 



24 Andreas H. Pawlik & Joop Schaye 






- 











Figure 16. Test 3: Effect of angular resolution. The same slice 
as shown in Fig. 1151 Contours show a neutral fraction of r] = 0.5 
(left-hand column) and r\ = 0.01 (right-hand column) at time 
t = 0.05 Myr (top row) and t = 0.2 Myr (bottom row). Green, red 
and blue lines correspond to the low (N c = 8), fiducial (N c = 32) 
and high (N c = 128) angular resolution simulations, resp. The 
numbers of ViPs per SPH particle at the end of the simulations 
are 0.05, 1 and 2.8 for the low, fiducial and high angular res- 
olution simulations, resp. Note that the fiducial simulation is al- 
ready converged, even though its angular resolution 7V C = 32 
corresponds to a relatively large cone opening angle of u> 41 
degrees. 



Figure 17. Test 3: Effect of resampling. The same slice as shown 
in Fig. 1151 Contours show a neutral fraction of r\ = 0.5 (left- 
hand column) and r] = 0.01 (right-hand column) at time t = 
0.05 Myr (top row) and t = 0.2 Myr (bottom row). The red 
contours correspond to the fiducial angular resolution simulation 
and are identical to the red contours in Figs. IT5land ll6l The blue 
contours correspond to a simulation identical to the simulation 
employing the fiducial angular resolution, except for the fact that 
in this simulation we periodically (every 10th radiative transfer 
time step) resampled the density field to suppress the particle 
noise. Note that resampling does not visibly decrease the effective 
angular resolution. 



torn row) for the fiducial angular resolution. The ionization 
front (neutral fraction r\ = 0.5) is delayed by the dense fila- 
ments, leading to the characteristic "butterfly" shapes of the 
ionized regions. For comparison, we also show the results ob- 
tained with two othe r codes, the ray-tracing scheme C 2 -RAY 
jMellema et alj|2006l; green contours) and the Monte Carlo 
code crash (jMaselli. Ferrara. fc Ciardil I200&1 ; ICiardi et"al] 
l200ll ; blue contours) , as published in the cosmological rad ia- 
tive transfer code comparison project (|lliev et al.l l2006 s R 
Both C 2 -RAY and CRASH are mesh codes, working directly 
on the uniform m esh input density field provided by the 
PM+TVD code of iRvu et all (|l993h - 

The agreement between the results of traphic and C 2 - 
RAY resp. CRASH is very good, not only for the ionization 
front, but also for smaller neutral fractions. The contours 
from traphic are slightly noisier than those from C 2 -RAY, 
which is expected since in addition to the particle noise 
affecting the radiative transfer, the Monte-Carlo sampling 
noise imprinted on the density field affects our simulations, 
particularly in the under-sampled low density regions, as 



15 The perfor mance of two more codes was reported in 
llliev et all (120061): FTTE l lRazoumov fe CardaHl2005T> and SIMPLEX 
llRitzerveld fc IckjlioOrl . For clarity and since they are very sim- 
ilar to the results obtained with C 2 -RAY and CRASH, we do not 
show the results obtained with FTTE. We do not include the re- 
sults of SIMPLEX in our comparison, since they differ considerably 
from those obtained with all other codes. 



already noted earlier. The noise level is, however, substan- 
tially lower than one would anticipate based on the tests 
presented in Sections 15.3.11 and 15.3.21 The most likely ex- 
planation for this welcome surprise is that the presence of 
multiple ionizing sources leads to a regularization in the dis- 
tribution of the neutral fraction. Numerical noise arising 
from the representation of the continuous density field by 
a discrete set of particles is therefore reduced. Differences 
between our results and those of C 2 -RAY resp. CRASH also 
arise through the different treatment of the photon spec- 
trum. Since the photo-ionization cross-section depends on 
frequency (Eq. I15[) . the thickness of finite ionization fronts 
(e.g. defined as 0.9 < 77 < 0.1) and hence the position of 
the particular contour r\ — 0.5 will in part be determined 
by the details of the numerical implementation of the multi- 
frequency transport. 

In Fig.[T5]we show the neutral fraction contours 77 = 0.5 
(left-hand side) and 77 = 0.01 (right-hand side) at times 
t = 0.05 Myr (top row) and t = 0.2 Myr (bottom row) for 
the low (green contours) and the high (blue contours) an- 
gular resolution simulations. For comparison, the contours 
obtained in the fiducial simulation are also shown (red) . We 
note that the high angular resolution simulation yields neu- 
tral fraction contours that are almost identical to those ob- 
tained in our fiducial simulation, indicating numerical con- 
vergence. The low angular resolution simulation, although 
still in good agreement with the high angular resolution 
simulation, fails to properly reproduce the expected neutral 
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fraction contours when scrutinised in detail. In the low angu- 
lar resolution simulation, neutral fraction contours are often 
slightly advanced instead of delayed by the dense filaments. 
Although the effect is small, it becomes apparent when the 
contours are compared to the corresponding contours of the 
high angular resolution simulation. 

The last observation agrees with the discussion of 
anisotropies in particle-to-neighbour transport schemes 
sketched in Section 14.2.11 and presented in detail in Ap- 
pendix [X] In Appendix |A1 we demonstrate that when pho- 
tons are transported from a given particle to its neighbours, 
the net transport direction is generally strongly correlated 
with the direction towards the centre of mass of the neigh- 
bouring particles. As a result, the transport is partly gov- 
erned by the spatial distribution of the SPH particles. For 
cosmological simulations this implies that photons are pref- 
erentially transported along dense filaments, traphic prop- 
agates photons in cones to overcome this problem. The cones 
confine the photons to the solid angles they were emitted 
into, ensuring a correct transport of radiation on the scale 
of the chosen angular resolution. If the angular resolution is 
chosen too low to properly resolve the structures in the SPH 
density field, the transport is no longer independent of the 
geometry of the SPH simulation and artefacts may occur, as 
we see in Fig. [16] for the low angular resolution simulation. 
Even with an angular resolution as low as N c = 8, however, 
the artefacts are small. Nevertheless, it is clear that in order 
to properly solve the radiative transfer equation, the angular 
resolution must be chosen high enough to establish numer- 
ical convergence. Fig. [16] shows that an angular resolution 
of 7V C = 32, which corresponds to a relatively large opening 
angle of u) ~ 41 degrees (Eq. [9}, is already converged. The 
reason why a relatively poor formal angular resolution suf- 
fices, is, as already noted in the discussion of the sharpness 
of shadows thrown by opaque obstacles in Section ["5.3.2l that 
the photon transport with traphic is intrinsically adaptive 
in angle. 

In Fig. [17] we present results for a simulation that used 
the resampling technique presented in Section lB3l but which 
was otherwise identical to the fiducial simulation presented 
above. The neutral hydrogen densities of the SPH particles 
were not re-calculated according to the perturbed positions 
resulting from the resampling, but kept constant to avoid 
additional scatter in the density. The resampling leads to a 
reduction in the particle noise, which is most apparent for 
the r\ = 0.01 contours, since in case of the r\ = 0.5 the noise 
is already low, as noted above. Note that the resampling 
does not noticeable degrade the angular resolution. 

In Fig. [18] we show the evolution of the mean (over 
all particles i) ionized fraction, both volume-weighted, i.e. 
Xv = J2ihiXi/12i h i> and mass- weighted, i.e. x™ = 
2Zi m iXi/ X^i m i- Again, the results obtained with traphic 
are in excellent agreement with the results obtained with C 2 - 
RAY and CRASH. For the latter two we obtained the mean 
ionized fraction by averaging the ionized fraction reported in 
the cosmological radiative transfer code comparison project 
(|lliev et al.ll2006l ) over all grid cells i, i.e. Xv = Y^i Xi/J2i 1 
and Xm = J2iPiXi/T,iPi- 

The ratio of mass-weighted and volume- weighted mean 
ionized fractions is at early times slightly larger for the 
low angular resolution simulation than for the high angular 
resolution simulation, as can be seen in the bottom panel 
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Figure 18. Test 3: The volume- and mass-weighted mean ionized 
fractions, XV and Xm, resp., averaged over the whole simulation 
box as a function of time, for the low, fiducial (without and with 
resampling of the density) and high angular resolution simulation, 
as indicated in the legend. All results fall nearly on top of each 
other. Differences in Xm/xv are OI1 ly noticeable when xv is small. 
For comparison, we also show the results o btained with C 2 -RAY 
and CRASH as reported in llliev et al. (2006). 

of Fig. 1181 This is another manifestation of the fact that 
particle-to-neighbour transport is generally not independent 
of the geometry of the SPH simulation, resulting in photons 
being preferentially transported into high (particle) density 
regions. It once more underlines the importance of the con- 
cept of emission and transmission cones (with sufficiently 
small solid angles) which traphic uses to accomplish the 
transport of radiation independently of the spatial distribu- 
tion of the SPH particles. 

Finally, we show that for the present radiative transfer 
problem simulations solving the time-dependent radiative 
transfer equation give a significantly different result than the 
simulations solving the time-independent radiative transfer 
equation discussed above. We carried out a simulation using 
N c = 32, N ng b — 32 and additionally employed the photon 
packet clocks to limit the distance over which photon packets 
can propagate during each time step. The size of the radia- 
tive transfer time step was set to At r = 10~ 3 Myr. This 
time step is a factor 10 times larger than the time step used 
for solving the time-independent radiative transfer equation 
in the simulations presented above, which required a smaller 
time step because of our particular treatment of the time- 
independent radiative transfer equation (see Section [5.2.3p . 
The locations of the ionization fronts (i.e. r\ — 0.5) ob- 
tained in this simulation are shown in Fig. 1191 fblue curves), 
together with those obtained in the corresponding simula- 
tion solving the time-independent radiative transfer equa- 
tion (red curves, cp. the left-hand panels of Fig. [15}, at four 
different simulation times t = 0.05, 0.1, 0.2 and 0.3 Myr. 

It is clear from Fig. [19] that the simulation solv- 
ing the time-independent radiative transfer equation pro- 
duces ionized spheres that are unphysically large. This is 
due to the well-known fact ( see, e.g., the discussion in 
I Abel. Norman, fc Madaulll999h that ionization fronts may 
propagate at speeds larger than the speed of light, if the 
time- independent radiative transfer equation is solved. The 
difference between the two simulations is larger at early 
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Figure 19. Test 3: Slice through the density field at z = Lt ox /2. Contours show a neutral fraction of r\ = 0.5 at times t = 0.05,0.1,0.2 
and 0.3 Myr (from left to right). Red contours show the results of our fiducial (7V C = 32, N ng h = 32) simulation (at times 0.05 and 
0.2 Myr these are identical to the red contours shown in the top-left resp. bottom-left panel of Fig. 1151 1. Blue contours show the result of 
a simulation employing the same resolution [N c = 32, N ng b = 32), but using the clock of the photon packets to solve the time-dependent 
radiative transfer equation. Note that in the simulation solving the time-independent radiative transfer equation the ionized regions are 
too large, since in this simulation the ionization fronts are initially propagating at speeds larger than the speed of light. 



times than at late times, which is expected, since in equi- 
librium, i.e. t r — > ao, the results of both simulations must 
agree. 

In summary, in this section we studied the propagation 
of ionization fronts around multiple sources in a static cos- 
mological density field. We demonstrated the importance 
of the concept of cones which underlies the photon trans- 
port in traphic. Without the confinement by transmission 
cones of sufficiently small solid angle, particle-to-neighbour 
transport is governed in part by the spatial distribution of 
the particles, resulting in the preferential transport of pho- 
tons into high (particle) density regions. Thanks to the fact 
that traphic is adaptive in angle, a relatively modest for- 
mal angular resolution of iV c = 32 is already sufficient to 
obtain a converged solution. Since the setup of our simula- 
tions followed the description of the corresponding test in 
the cosmological radiative transfer code comparison project 
ijlliev et al.ll2006t ). we were able to benchmark our radiative 
transfer scheme by comparing wi th the results obtain ed by 



the ray-tracing scheme C 2 -RAY (Mcllcma ct al. 200 tj) and 



the Monte Carlo code crash l|Maselli. Ferrara. fc Ciardil 
120031 ; I Ciardi et al.ll200lh . We found excellent agreement in 
the positions of neutral fraction contours as well as the mass 
and volume- weighted mean ionized fractions. 

We have furthermore seen that for the test problem 
presented in this section, simulations solving the time- 
independent radiative transfer equation lead to ionized re- 
gions that are unphysically large during their early evolu- 
tion. This illustrates the importance of correctly accounting 
for the finite speed of light when performing radiative trans- 
fer simulations to study the morphology of ionized regions 
that are strongly out of equilibrium. 



6 CONCLUSIONS 

In this paper we have presented traphic (TRAnsport of 
PHotons In Cones), a novel radiative transfer scheme for 
SPH simulations, traphic has been designed for use in 
radiation-hydrodynamical simulations exhibiting a wide dy- 
namic range of physical length scales and containing a large 
number of light sources, as for instance encountered in cos- 
mological simulations. The requirements for this are high: 



Due to the limits on the computational power available to- 
day, radiative transfer simulations at the resolution of state- 
of-the-art hydrodynamical simulations need to be adaptive 
both in space and in angle and parallel on distributed mem- 
ory machines. 

traphic meets these requirements. The radiative trans- 
fer equation is solved by transporting photons in an explic- 
itly photon-conserving manner directly on the discrete set 
of SPH particles. Photons are transported globally by prop- 
agating them locally, between a particle and its neighbours. 
For a fixed number of SPH particles, the number of neigh- 
bours employed for the transport, the parameter N ng b, sets 
the adaptive spatial resolution. The conceptual similarity 
between the photon transport and the calculation of parti- 
cle properties in SPH simulations allows a straightforward 
numerical implementation of traphic in SPH simulations 
for application on distributed memory machines. 

In traphic photon packets are transported inside 
cones. Because source particles emit into a set of cones that 
tessellates the simulation box, the angular dependence of 
the emission can be modelled independently of the spatial 
distribution of the neighbouring SPH particles. The number 
of cones, N c , is the second parameter employed in traphic 
and determines the formal angular resolution of the radia- 
tive transfer simulation. After their emission, the photon 
packets are transported downstream from SPH particle to 
SPH particle. They are confined to the solid angle they 
were originally emitted into by regular cones of solid an- 
gle 47r/jV c . Mimicking the ray-splitting technique used in 
ray-tracing schemes, the photon transport is adaptive in an- 
gle. The effective angular resolution is therefore much higher 
than the formal one. Radiative transfer simulations employ- 
ing traphic can thus be performed using only a relatively 
small number of cones without loss of accuracy. 

The propagation of photons inside cones that do not 
contain a neighbour requires the introduction of so-called 
virtual particles (ViPs). It is the concept of cones combined 
with ViPs that enables the directed transport of photons on 
the unstructured grid defined by the discrete set of irregu- 
larly distributed SPH particles. In addition to the number of 
cones N c and the number of SPH neighbours N ng b, it is the 
ratio N c /N ng b that directly influences the performance of 
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traphic. It controls the amount of noise introduced by the 
representation of the underlying continuum physics with a 
discrete set of particles. This particle noise is small for both 
N c < N ng b and N c > N ng b due to the large number of 
neighbours per cone and the large number of ViPs, resp. 
For the choice of parameters N c — N ng b the particle noise 
can be substantial. It can, however, be efficiently suppressed 
by employing a density field resampling strategy. 

A practical problem often encountered when solving the 
radiative transfer equation is the linear scaling of the com- 
putational effort with the number of sources. Such a scaling 
limits the application of most of the existing radiative trans- 
fer schemes to simulations containing only a few sources. 
In traphic this problem is circumvented by introducing a 
source merging procedure. Photon packets emitted by dif- 
ferent sources are merged in accordance with the angular 
resolution employed, such that at any point in space at most 
N c photon packets need to be propagated. This renders the 
photon transport effectively independent of the number of 
sources in the simulation and makes it feasible to include 
a diffuse radiation component emitted by the SPH parti- 
cles. In the most extreme case of radiation from an arbitrary 
number of sources completely filling the simulation box, the 
computational effort required by traphic merely scales as 
the product of spatial and angular resolution, N x N ng b x N c , 
where JV is the total number of particles (SPH + stars). 

Applications of traphic include the solution of both 
the time-independent and time-dependent radiative transfer 
equation in large-scale simulations of cosmological reioniza- 
tion. Here we have specified traphic for the transport of 
mono-chromatic hydrogen-ionizing radiation and described 
its implementa tion in the parallel SPH code GADGET-2 
|Springell l2005). We showed that traphic is able to accu- 
rately reproduce the expected growth of the ionized sphere 
around a single point source in a homogeneous medium and 
to cast sharp shadows behind opaque obstacles. Further- 
more, we tested our scheme in a physical setting of complex 
geometry: The growth of ionized regions around multiple 
point sources in a cosmological density field. De tailed com- 
paris ons with the results obtained by other codes (|lliev et al.l 
2006) for identical problems show excellent agreement. 

In this paper we have limited ourselves to the descrip- 
tion of the radiative transfer scheme traphic and its imple- 
mentation for use on static density fields. However, traphic 
can by design readily be coupled to the thermo- and hydro- 
dynamic evolution of matter in SPH simulations. We will 
report on such a coupling and its extension to the transport 
of multi-frequency radiation in future work. 



ACKNOWLEDGMENTS 

We thank Claudio Dalla Vecchia for stimulating discussions, 
enduring encouragement as well as technical help and Huub 
Rottgering for his valuable advice. We are grateful to Craig 
Booth and Benedetta Ciardi for a thorough reading of the 
draft. A.H.P. thanks the University of Potsdam and the Max 
Planck Institut fur Astrophysik for their hospitality. Some 
of the simulations presented here were run on the Cosmol- 
ogy Machine at the Institute for Computational Cosmology 
in Durham as part of the Virgo Consortium research pro- 
gramme and on Stella, the LOFAR BlueGene/L system in 



Groningen. This work was supported by Marie Curie Excel- 
lence Grant MEXT-CT-2004-0I4II2. 



REFERENCES 

Abel T., Norman M. L., Madau P., 1999, ApJ, 523, 66 
Abel T., Wandelt B. D., 2002, MNRAS, 330, L53 
Alvarez M. A., Bromm V., Shapiro P. R., 2006, ApJ, 639, 
621 

Barkana R., Loeb A., 2004, ApJ, 609, 474 
Brookshaw L., 1994, MmSAI, 65, 1033 
Cen R., 2002, ApJS, 141, 211 

Ciardi B., Ferrara A., Marri S., Raimondo G, 2001, MN- 
RAS, 324, 381 

Croft R. A. C, Altay G, 2007, arXiv, 709 JarXiv:0709.2362l 
Dale J. E., Ercolano B., Clarke C. J., 2007, MNRAS, 1056 
Dale J. E., Bonnell I. A., Whitworth A. P., 2007, MNRAS, 
375, 1291 

Dopita M. A., Sutherland R. S., 2003, Astrophysics of the 

diffuse universe, Springer 
Fryer C. L., Rockefeller G, Warren M. S., 2006, ApJ, 643, 

292 

Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375 
Gingold R. A., Monaghan J. J., 1978, MNRAS, 184, 481 
Gingold R. A., Monaghan J. J., 1982, J. Comput. Phys., 
46, 429 

Gnedin N. Y., Abel T., 2001, New A, 6, 437 
Goldstein H., 1980, Classical Mechanics, Addison- Wesley 
Publ. Co 

Gritschneder M., Naab T., Heitsch F., Burkert A., 2007, 
IAUS, 237, 246 

Herant M., Benz W., Hix W. R., Fryer C. L., Colgate S. A., 
1994, ApJ, 435, 339 

Hernquist L., Katz N., 1989, ApJS, 70, 419 

Hockney R. W., Eastwood J. W., 1988, Computer Simula- 
tions Using Particles, Taylor & Francis 

Iliev I. T., Mellema G, Pen U.-L., Merz H., Shapiro P. R., 
Alvarez M. A., 2006, MNRAS, 369, 1625 

Iliev I. T., et al., 2006, MNRAS, 371, 1057 

Johnson J. L., Greif T. H., Bromm V., 2007, ApJ, 665, 85 

Kessel-Deynet O., Burkert A., 2000, MNRAS, 315, 713 

Kohler K., Gnedin N. Y., Hamilton A. J. S., 2007, ApJ, 
657, 15 

Lucy L. B., 1977, AJ, 82, 1013 

Maselli A., Ferrara A., Ciardi B., 2003, MNRAS, 345, 379 
Mayer L., Lufkin G, Quinn T., Wadsley J., 2007, ApJ, 661, 
L77 

Mellema G, Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, 
New A, 11, 374 

Mihalas D., Weibel Mihalas B., 1984, Foundations of radi- 
ation hydrodynamics, Oxford University Press, New York 
Miles R. E., 1965, Biometrika, Vol. 52, No. 3/4, pp. 636 
Monaghan J. J., 1992, ARA&A, 30, 543 
Monaghan J. J., Rep. Prog. Phys., 68, 1703 (2005) 
Okabe A., Boots B., Sugihara K., Chiu S. N., 2000, Spatial 

Tessellations, Second edition, Wiley 
Osterbrock D. E., 1989, Astrophysics of gaseous nebulae 

and active galactic nuclei, Palgrave Macmillan 
Oxley S., Woolfson M. M., 2003, MNRAS, 343, 900 
Paschos P., Norman M. L., Bordner J. O., Harkness R., 
2007, arXiv, 711 JarXiv:0711~1 904 



28 Andreas H. Pawlik & Joop Schaye 



Press W. H., Teukolsky S. A., Vetterling W. T., Flannery 
B. P., 1992, Numerical recipes in C. The art of scientific 
computing, Cambridge University Press, 2nd ed. 
Price D., 2005, astro, [arXiv:astro-ph/0507472| 
Razoumov A. O., Cardall C. Y., 2005, MNRAS, 362, 1413 
Razoumov A. O., Sommer-Larsen J., 2006, ApJ, 651, L89 
Ritzerveld J., Icke V., 2006, PhRvE, 74, 026704 
Ryu D., Ostriker J. P., Kang H., Cen R., 1993, ApJ, 414, 
1 

Semelin B., Combes F., Baek S., 2007, arXiv, 707, 
larXiv:0707.2483l 

Soneira R. M., Peebles P. J. E., 1978, AJ, 83, 845 
Springel V., Hernquist L., 2002, MNRAS, 333, 649 
Springel V., 2005, MNRAS, 364, 1105 
Stamatellos D., Whitworth A. P., 2005, A&A, 439, 153 
Stewart G. W., 1980, SIAM Journal on Numerical Analysis, 

Vol. 17, No. 3, 403 
Susa H., 2006, PASJ, 58, 445 

Susa H., Umemura M., 2006, ApJ, 645, L93 

Trac H., Cen R., 2006, astro, arXiv:astro-ph/0612406"1 
Viau S., Bastien P., Cha S.-H., 2006, ApJ, 639, 559 
Whalen D., Norman M. L., 2007, arXiv, 708, 
larXiv:0708.2444l 

Wise J. H., Abel T., 2007, arXiv, 710, larXiv:0710.3160l 
White S. D. M., 1996, Cosmology and large scale structure, 
Proceedings of the "Les Houches Ecole d'Ete de Physique 
Theorique" (Les Houches Summer School), p. 349, Else- 
vier Scientific Publishing Company, Amsterdam 
Whitehouse S. C, Bate M. R., 2006, MNRAS, 367, 32 
Yoshida N., Oh S. P., Kitayama T., Hernquist L., 2007, 
ApJ, 663, 687 



APPENDIX A: THE ANISOTROPY OF 
PARTICLE-TO-NEIGHBOR TRANSPORT 

In Section ^. 2. ll we introduced an emission cone tessellation 
to accomplish the emission of photons by a source particle to 
its neighboring gas particles. Then we also mentioned that 
the cones are required to achieve an emission in agreement 
with the angular dependence of the emissivity of the source, 
independently of the spatial distribution of the neighbors. 
In this appendix we explain why this is the case. In par- 
ticular, assuming that particles have equal mass, we show 
that when transporting photons (or more generally, any ex- 
tensive quantity) in particle simulations from a particle to 
its neighbors, the net transport direction is generally corre- 
lated with the direction towards the centre of mass of the 
neighbouring particles. As a result, the transport is partly 
governed by the geometry of the simulation, in addition to 
the intrinsic emissivity of the source. 

Consider a particle located at the origin O of a 3- 
dimensional coordinate system (the emitter frame) . The par- 
ticle emits photons to its N ng b nearest neighbours, residing 
in the sphere of radius h centred on the emitting particle. 
Throughout this section we assume that all particles have 
equal mass, m = 1. The emission is performed by transfer- 
ring the fraction Wi / ^^T-f b Wj of the total number of pho- 
tons to neighbour i (i = 1, 2, N ng t), where the Wi ^ are 
weights to be specified. The emission properties will depend 



both on the weights and on the spatial distribution of the 
neighbours. 

To study the angular dependence of the emission pro- 
cess, we introduce the vector sum s = J^. WiUi/ £\ Wi over 
all unit vectors = r»/|ri|, where r; is the position of 
neighbour i. This vector sum can be interpreted as the net 
direction of the emission. Let us denote the angle enclosed 
by s and the vector to the centre of mass r cm of the neigh- 
bours, r cm = Vi/ N ng b, with a. We ask for the probability 
density function (pdf) p(cosa) of the cosine of a, where 

s-r cm >l "Mi. ■ }] r, 

COSQ=,^ r = t= rp= (Al) 

|s||r cm | |Y\ u^Ui 1 1 2Zj r A 

The meaning of the pdf can be understood by look- 
ing at its extremes. If p(cosa) is strongly peaked around 
cos a ~ 1, the emission has a net direction biased towards 
the centre of mass direction, and photons will be preferen- 
tially transported into the high (particle) density regions. 
On the other hand, if p(cos a) is strongly peaked around 
cos a w — 1, photons will be preferentially transported into 
the direction opposite to the centre of mass. Finally, a flat 
pdf p(cosq) = 1/2 indicates that the net emission direction 
and the direction towards the centre of mass are uncorre- 
cted. 

As an illustration, assume that all neighbours have iden- 
tical weights Wi = 1 and that their positions r; are drawn 
from a probability distribution that uniformly samples the 
surface of the sphere of radius h surrounding the emitter at 
O. To obtain p(cosa), we note that our assumptions imply 
that the centre of mass vector r cm and the vector sum s point 
in the same direction, so that p(cosa) = 25d(cosq — 1), 
where 5d{x) is the Dirac S -function. Hence, there is a per- 
fect correlation between s and r cm . 

Another example is given by assuming that the particles 
of unit mass are distributed randomly within the sphere of 
radius h around the emitter. The resulting pdf p(cosa) as 
obtained through a Monte Carlo simulation is shown in its 
cumulative form F(cosa) = f dx p(x) in Fig. IA1I for 

two choices of the number of neighbours, N ng b = 16 (thin 
dotted curve) and N ng b = 64 (thick dotted curve; falling 
nearly on top of the thin dotted curve). Clearly, there is a 
very strong correlation between the net direction of emission 
and the direction towards the centre of mass of the neighbour 
distribution. 

For the important case of isotropic sources of photons 
the net emission direction should not be correlated with the 
direction towards the centre of mass, i.e. p(cosa) = 1/2, 
independent of cos a, and thus F(cosq) = (1 + cosa)/2. 
From here on, we will refer to this case as isotropic emission 
of radiation. We emphasize the statistical character of this 
statement. An individual particle may still emit anisotrop- 
ically, which could be revealed by studying the amplitude 
|s|. For emission that is isotropic for individual particles, 
|s| = 0. Consider the neighbour distribution of the last ex- 
ample, for which we showed the emission pattern to be far 
from isotropic. Can the emission be isotropized by tuning 
the weights? 

For instance, consider solid angle weighting, defined as 
follows. We project all neighbours radially onto the unit 
sphere. To each neighbour we attach a weight proportional 
to its area of influence, which we define by the collection of 
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all points on the unit sphere which are closer to the pro- 
jected neighbour for which we calculate the weight than to 
all other projected neighbours. Hereby, the distance of two 
points on the sphere is given by the arclength of the con- 
necting geodesic (great circle). This definition of the area of 
influence results in a specific tessellation of the surface of the 
spher e, the so-called Voronoi tessellation (e.g. lOkabe et al.l 
|200J ). and provides us with well-defined non-overlapping 
solid angles. 

The resulting cumulative distribution _F(cos a) is shown 
m Fig. |AT] (dashed curves) . Although the weighted emission 
has weakened the correlation between s and r cm as com- 
pared to the case of equal weights, the net emission direction 
is still markedly peaked towards the direction to the centre 
of mass. We note that weighting by solid angle is the most 
natural choice and its inability to isotropize the transport 
might be surprising. Alternatively, one could introduce ad 
hoc weights w% = vol, where r > 1 is some exponent and 
the Wi are solid angle weights, to re-shape the distribution 
p(cosq) to become less and less peaked around cos a ~ 1. 

However, the reason why the solid angle weights fail to 
isotropize the transport is not that the weights are inappro- 
priate, but that the directions to the neighbours are fixed 
and thus cannot be freely chosen along with the weights, as 
mentioned in Section [4.2.11 This strongly limits one's ability 
to influence the shape of p(cosa). For illustration, consider 
an emitter having all its neighbours located in one hemi- 
sphere. The net emission direction s will then necessarily lie 
in that hemisphere, too, regardless of the chosen weights. 
For a statistical ensemble of emitters, s will then necessarily 
be correlated with the direction towards the centre of mass. 

This also explains the dependence of p(cosa) on the 
number of neighbours. The larger N ng b, the greater the 
probability of finding a neighbour in any chosen solid an- 
gle. Through this increase in directional sampling, p(cosct) 
can be expected to approach the uncorrelated case. This 
behaviour can indeed be observed in Fig. IA1I 

We also compute the statistic p(coscv) for the case of 
neighbours clustered in space, which is more typical for cos- 
mological simulations of structure formation. We did this 
both using a toy clustering model (|Soneira fc Peeblesll 19781 ) 
and for a cosmological density field obtained from a ACDM 
hydrodynamical simulation (the same field that is used in 
Section r5.3.3[) . For the weighting schemes we investigated the 
dependence of the shape of p(cos a) on the chosen weights is 
qualitatively very similar to that in the case of the random 
distribution of neighbours discussed earlier. 

We emphasize that our discussion is not limited to the 
transport of photons. Indeed, particle-to-neighbour trans- 
port is for example routinely employed in SPH simulations 
of galaxy formation, where a star particle distributes its met- 
als and energy amongst its neighbours. Different recipes for 
distributing metals and energy can be found in the litera- 
ture. Most of these apply a weighting related to the SPH 
formalism. As an example, to each neighbour i we associate 
the weight 



(A2) 



1.0 
0.8 

«T 0.6 

ID 

o 
o 

~" 0.4 

0.2 
0.0 



SPH NG64 
SPH NG16 
Solidongle NG64 
Solidangle NG16 
Equal NG64 
Equal NG16 




1.0 



-0.5 0.0 

cos a 



0.5 



.0 



Here, h is the radius of the neighbourhood of the emitting 
particle (at O), Ti = |r»| and W is the spline kernel Eq. [4] 



Figure Al. Monte Carlo simulation of particle-to- neighbour 
transport. We show the cumulative probability distribution of 
cos a, where a is the angle between the net emission direction s 
and the direction towards the centre of mass r cm of the neigh- 
bours. The curves correspond to the following weighting schemes, 
with the number of neighbours indicated in brackets (from top to 
bottom): SPH (64), SPH (16), Solidangle (64), Solidangle (16), 
Equal (64), Equal (16). The last two curves are nearly on top of 
each other. 



Performing the transport on the same random neighbour 
distribution as before, we obtain the distribution p(cosa). 
As can be seen in Fig. IA1I (solid curves) , the direction of net 
transport and the direction to the centre of mass are found to 
be almost uncorrelated. We note that this result, the reason 
of which at the moment is not clear to us, is only a statisti- 
cal statement. We calculated |s| and found that the emission 
of individual particles is still not close to isotropic. We ar- 
rive at qualitatively similar results for different definitions 
of the SPH weights, both for random and clustered neigh- 
bour distributions. A more quantitative statement, however, 
must depend on both the specific properties of the spatial 
distribution of the neighbours, and the precise form of the 
adopted SPH weights. 

In summary, in this appendix we have shown that when 
transporting a quantity from a particle to its neighbours, 
the net transport direction is generally governed by the 
spatial distribution of the neighbours, in addition to the 
the intrinsic properties of the emitting particle. We empha- 
size that our observations apply to the particle-to-neighbour 
transport of arbitrary quantities in any particle simulation. 
We would like to remind the reader that in Section [3] we 
presented a specifically designed particle-based transport 
scheme, traphic, that overcomes the problems of particle- 
to-neighbour based transport discussed in this appendix, not 
just statistically, but also on a particle by particle basis. We 
achieved this by employing cones to confine photons to the 
solid angle into which they are emitted, making the trans- 
port independent of the geometry of the SPH simulation, on 
the scale of the chosen angular resolution. 
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APPENDIX B: NUMERICAL 
IMPLEMENTATION 

Bl Cones 

In this section we describe the numerical implementation of 
the cones employed for the emission and reception of pho- 
ton packets in traphic. For the emission (Section 14.2. 
each source particle divides its neighbourhood using a set 
of tessellating emission cones. The same tessellation is also 
employed for the reception of photon packets by gas par- 
ticles (Section I4.2.3|) . In the following, we employ spherical 
coordinates (r, <f>, 0), which are related to the Cartesian com- 
ponents (r x ,r y ,r z ) of an arbitrary vector r through 

(Bl) 
(B2) 
(B3) 

In our implementation, the emission (reception) cones sam- 
ple the volume around each particle isotropically. Since 
the surface element of the unit sphere is given by ds — 
d(cos9)d(j>, this is achieved by distributing the cone bound- 
arie^f] uniformly (i.e. on a regular lattice with indices i, j) in 
(cos 0, (f>). Thus, the boundaries of cone are described 

by the four arcs of constant 
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Correspondingly, we define the emission (reception) cone 
axes by 
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Note that each of the N c — N$ X No emission (recep- 
tion) cones has the same solid angle fl = 4n/N c , as can 
be seen from integrating over the surface element of the 
unit sphere within the boundaries (Eqs. IB4IIB7[) of a sin- 
gle cone. We also implemente d the tessellation used in 
lAbel. Norman, fc MadavJ (|l999l ). which leads to cones that 
are more similar in shape. We could not find any systematic 
differences in the test problems described in Sections 15.3.11 
- 15.3.31 using either of the two types of tessellations. This 
is not surprising because any artefacts due to the shape of 
the cones will be suppressed by the random rotations of the 
emission (reception) cones we perform before each emission 
(reception), as we will describe in the following section. 



tessellation has a random orientation. The primary motiva- 
tion for randomly rotating cones is to increase the angular 
sampling. Furthermore, randomly rotating the emission and 
reception cones leads to a reduction of artefacts arising from 
the particular definition we employ to construct the cone tes- 
sellation, as noted in the last section. Here we describe our 
numerical implementation of the random rotations. 

We can think of the set of cones that comprises a partic- 
ular cone tessellation as a rigid body, to which we can attach 
a local Cartesian coordinate system with axes {e' x , e' y , e' z } . 
The orientation of this coordinate system with respect to the 
canonical Cartesian coordinate system, e.g. the simulation 
box axes {e x ,e y ,e z }, is fully described by three variables, 
the Eulerian angles (e.g. lGoldstein|[l980l ). 

Eulerian angles are defined as the three successive an- 
gles of rotations that map the axes {e x ,e y ,e z } onto the 
axes {e' x ,e y ,e'z}. There exist several conventions. In the 
zxz convention we employ here, the initial system of axes 
{e x ,e y ,e z } is first rotated counter-clockwise by an angle cfi 
around the z-axis, with the resulting coordinate system la- 
belled {e^, e v , }. Second, the coordinate system {e^,e v , e^} 
is rotated by an angle 6 counter-clockwise about the £-axis, 
leaving the new coordinate system {e'^,e' v ,e^}. The third 
and last rotation is carried out counter-clockwise by an an- 
gle tp around the ("'-axis, giving the desired {e' x , e' y , e' z } co- 
ordinate system. 

To obtain random Eulerian angles, we note that the 
invariant measure \x (the "volume element") on SO (3), the 
group of proper rotations i n R 3 , in th e zxz Eulerian angle 
parametrization reads (e.g. lMileall965l ). 
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Random Eulerian angles are therefore obtained by draw- 
ing random variables U\,U2, 113 from a uniform distribution 
on t he interval [0, 1] a nd applying the usual transformation 
(cp. iPress et al.ll 19921 ) . 



cj> = 2tyu\ 

9 = arccos(l — 2u2) 

if) = 2TTU:j. 
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We implement random rotations using rotation matrices, 
which are obtained from the random Eulerian angles. The 
matrix elements of a matrix representing a rotation associ- 
ated with a given set of Eulerian angles c an be readily ca l- 
culated. For details, we refer the reader to lGoldst ein ( 1980). 

In principle, random rotation matrices can be obtained 
witho ut drawing random Eulerian angles (e.g. IStewartj 
1980). However, we find that it is faster to generate ran- 
dom Eulerian angles, and then calculate the corresponding 
rotation matrices. Moreover, storing the three Eulerian an- 
gles instead of the nine rotation matrix elements reduces the 
memory cost. 



B2 Random rotations 

Recall from Sections l4.2.2l and l4.2.3l that we apply a random 
rotation to each cone tessellation. Consequently, each cone 

16 Strictly speaking, one should distribute the cone axes uni- 
formly in (cos but this implies asymmetric cones. 



B3 Reduction of particle noise 

For some of the simulations presented in Sections 15.3.11 - 
15.3.31 we employ the resampling approach described in Sec- 
tion 33] t° reduce numerical artefacts arising from the rep- 
resentation of the continuous density field with a discrete 
set of particles. The resampling requires the sampling of the 
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SPH kernel used in GADGET-2, which is the spline given by 
Eq. [4] We approximate this kernel by a Gaussian, 

W ° {r)= (2^W eXp( - r2/2<T2) - (B14) 

We find that with a standard deviation of a ~ 0.29 x h, the 
Gaussian provides a rea sonable fit to the spline. A similar re - 
lation was employed in lAlvarez. Bromm. fc Shapiro! (2006). 
When we resample the SPH density field, all SPH particles 
are redistributed by randomly displacing them from the po- 
sitions given by the hydrodynamical simulation, with the 
probability to find a given particle displaced by the distance 
r given by Eq. IB14I The (rare) displacements larger than h 
are discarded; in this case the original particle positions as 
determined by the hydrodynamical simulation are used. 
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